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Abstract This paper is concerned with the design, analysis and implementation of preconditioning con- 
cepts for spectral DG discretizations of elliptic boundary value problems. The far term goal is to obtain 
robust solvers for the "fully flexible" case. By this we mean Discontinuous Galerkin schemes on locally re- 
fined quadrilateral or hexahedral partitions with hanging nodes and variable polynomial degrees that could, 
in principle, be arbitrarily large only subject to some weak grading constraints. In this paper, as a first 
step, we focus on varying arbitrarily large degrees while keeping the mesh geometrically conforming since 
this will be seen to exhibit already some essential obstructions. The conceptual foundation of the envisaged 
preconditioners is the auxiliary space method, or in fact, an iterated variant of it. The main conceptual 
pillars that will be shown in this framework to yield "optimal" preconditioners are Legendre-Gaufi-Lobatto 
grids in connection with certain associated anisotropic nested dyadic grids. Here "optimal" means that 
the preconditioned systems exhibit uniformly bounded condition numbers. Moreover, the preconditioners 
have a modular form that facilitates somewhat simplified partial realizations at the expense of a moderate 
loss of efficiency. Our analysis is complemented by careful quantitative experimental studies of the main 
components. 

AMS subject classification 65N35, 65N55, 65N30, 65N22, 65F10, 65F08 

Keywords Discontinuous Galerkin discretization for elliptic problems, interior penalty method, locally 
refined meshes, variable polynomial degrees, auxiliary space method, Gaufi-Lobatto grids, dyadic grids. 

I Introduction 

Attractive features of discontinuous Galerkin (DG) discretizations are on one hand their versatility 
regarding a variety of different problem types as well as on the other hand their flexibility regarding 
local mesh refinement and even locally varying the order of the discretization. While initially the 
main focus has been on transport problems like hyperbolic conservation laws an increased attention 
has recently been paid to diffusion problems which naturally enter the picture in more complex 
applications like the compressible or incompressible Navier-Stokes equations. Well-posedness and 
stability issues are by now fairly well understood [H |3] although these studies refer primarily to 
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uniform fixed degrees and conforming domain partitions, see, e.g., |HLI12j for locally refined meshes 
and variable degrees. 

Of course, the option of choosing very high polynomial degrees is of particular interest for diffusion 
problems since under certain circumstances highly regular, even analytic solutions occur that are 
best treated by spectral discretizations |15j . The fact that in most practical situations a very high 
regularity is encountered only in part of the domain has led to the concept of /ip-discretizations which 
typically come about as conforming methods. However, the DG framework seems to be particularly 
well suited for accommodating strongly varying polynomial degrees (see, e.g., [241126)1 in combination 
with local mesh refinements admitting hanging nodes. 

While stability of such discretizations can be ensured, as will be seen below, essentially follow- 
ing standard lines, the issue of efficiently solving the resulting linear systems of equations seems to 
be much less clear. For quasi-uniform meshes and uniform (low) polynomial degrees a multigrid 
scheme, proposed in |18j . (covering actually convection-diffusion problems as well) does give rise 
to uniformly bounded condition numbers for diffusive problems provided that (i) the underlying 
hierarchy of meshes is quasi-uniform and (ii) the solution exhibits a certain (weak) regularity. This 
scheme has been extended in |19j to locally refined meshes showing a similar performance without a 
theoretical underpinning though. Domain decomposition preconditioners investigated in [HI2], give 
rise to only moderately growing condition numbers. A two- level scheme in the sense of the auxiliary 
space method (see, e.g., [3[22J[25]) nas been proposed in [T7] and shown to exhibit mesh-independent 
convergence again on quasi-uniform conforming meshes with a fixed uniform polynomial degree. In 
the framework of the auxiliary space method preconditioners providing uniformly bounded condi- 
tion numbers for locally refined meshes under weak grading constraints and variable but uniformly 
bounded polynomial degrees have been developed in [HJII2], but again the results are valid only for 
a fixed bound on the polynomial degrees with constants depending on that bound. 

The state of the art concerning preconditioners that are robust with respect to the polynomial 
degrees is somewhat heterogeneous. Meanwhile a fairly good understanding has been developed for 
pure spectral discretizations, i.e. for a single polynomial patch Q3J [231 E] ■ In the DG context the 
current results seem to draw primarily from domain decomposition concepts, see, e.g., |25j . To our 
knowledge they are subject to two type of constraints, namely to geometrically conforming meshes 
and to uniform poynomial degrees p. The best known results seem to offer bounds on condition 
numbers of only logarithmic growth in p. 

The objective of this paper is therefore to develop and analyze preconditioners for what we call 
Spectral- Element- DG discretizations, eventually exploiting the full flexibility potential of the DG 
concept, meaning arbitrarily large varying polynomial degrees and local refinements with hanging 
nodes, both only subject to mild grading constraints. The central theoretical challenge hings on 
the question how to design preconditioners that give rise to uniformly bounded condition numbers 
in the most flexible setting described above. We shall address this challenge in the context of the 
auxiliary space method [22J |27J [2H]- However, to keep the paper at an acceptable length, as a 
first step we confine the discussion here on geometrically conforming meshes, i.e., hanging nodes 
are avoided but varying polynomial degrees of arbitrary size are permitted. In fact, some essential 
obstructions are already encountered in this case for the following reason. First, on the basis of the 
experience in the context of spectral methods, in our opinion the most promising avenue towards 
fully robust preconditioners with respect to polynomial degrees is to employ concepts centering 
around Legendre-Gaufi-Lobatto grids (LGL-grids). As soon as the degrees on adjacent polynomial 
patches differ the corresponding LGL-grids do not match at the patch interfaces (the more so when 
dealing with hanging nodes) since LGL-grids are not nested. This turns out to have far reaching 
consequences in several respects, namely properly dealing with jump terms on the patch interfaces 
but also regarding the efficient solution of the patchwise problems in case of very large degrees. One 
of the remedies we propose is to employ auxiliary spaces based on certain anisotropic dyadic grids 
which are associated with the LGL-grids but which in addition are nested. Therefore, the concepts 
developed for the geometrically conforming case will be similarly relevant in the general case as well. 
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Aside from these principal theoretical issues the second main theme of the present studies is to 
explore the quantitative performance of the main building blocks entering such preconditioners. 

The layout of the paper is as follows. After formulating in Sect. [2] our model problem and the 
basic ingredients concerning Legendre-GauB-Lobatto (LGL)-grids, we recall in Sect. [3] the relevant 
conditions of the auxiliary space method. Sect. [4] offers a brief orientation concerning the subsequent 
development of a two-stage preconditioner. The first stage, using high-order conforming polynomials 
as auxiliary space is treated in Sect. [5] As a preparation of the subsequent second stage we discuss 
in Sect. [6] a major ingredient of our approach, namely certain hierarchies of nested dyadic grids 
associated with LGL-grids in a way that mutual interpolation between low order finite elements on 
these grids is stable in L 2 and H 1 . Based on these findings we develop and analyze in Sect. [7] stage 
II of our preconditioner where the new auxiliary spaces are built by low-order finite elements on 
dyadic grids. We conclude in Sect. [8] with detailed numerical studies of the various ingredients of 
the preconditioner to assess its quantitative performance in dependence of the various parameters 
entering the scheme. In particular, in Sect. |8.3TT] we discuss efficient ways of applying the "relaxation" 
operator defined by the auxiliary space method. 

Throughout the paper we shall employ the following notational convention. By a < b we mean 
that the quantity a can be bounded by fixed constant multiple of 6, independent on the parameters 
a and b may depend on. Likewise a ~ b means a < b and b < a. 



2 Formulation of SE-DG schemes 

In this section, we introduce the basic high-order DG discretization of a second-order model problem, 
as well as several variants of it. We start by recalling some essential concepts about Legendre methods 
and Legendre Gaufi-Lobatto grids on a single Cartesian element. 

2.1 (Hyper-) rectangles, polynomial spaces and LGL grids 

Let I — [a, b] be any bounded interval, of size H = b — a. For any fixed integer p > 1, let P p (7) be 
the space of the restrictions to / of the algebraic polynomials of degree < p on the real line. Let 
a = £ < • • • < < £j < ■ ■ ■ < £ p = b be the Legendre-Gaufi-Lobatto (LGL) quadrature nodes of 
order p in /, for which there exist weights Wo < • • • < Wj-i < Wj < • • • < w p such that 

p r 

]T t v(£j)w j = v(x)dx VweP 2p _i(I). (2.1) 

3=0 Jl 

We will denote by G P (I) the set of such nodes, which will be referred to as the LGL grid in /. 
Obviously, any v £ V p (I) is uniquely determined by its values on G P (I). We recall that nodes 
and weights are classically defined on the reference interval I — [—1,1], as uij, resp.. By affine 
transformation, one has £j = a + H(£j + l)/2 and Wj = (H/2)wj. We also recall that weights with 
index j close to or p satisfy Wj ~ Hp~ 2 (in particular, w = w p = H(p(p+ 1)) _1 ), whereas weights 
with index j close to p/2 satisfy Wj ~ Hp~ x . Yet, the variation in the order of magnitude is smooth, 
as made precise by the estimates 

1 < min < max < i ; (2.2) 

i<i<p Wj i<j<p Wj 

which hold uniformly in j, p and H (see, e.g., JEHUS])- A fundamental property for the sequel (see, 
e.g., [15 ) is that the bilinear form 



p 

{U, u)o,Z,p = U (Q V (& W 3 
3=0 



(2.3) 
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is an inner product in P p (7), which induces therein a norm uniformly equivalent to the L 2 -norm: 

IMkr < IHIo,/, P < V3|H| ,/ VveP p (I). (2.4) 

Furthermore, we shall make frequent use of the following "trace" inequality: for any v £ V p (I), one 
has 

Ke)|<^±i|M|o,/, (2.5) 



where e £ {a, 6} denotes any endpoint of /. This follows, by a scaling argument, from the expansion 
v(x) — ^2 P k = o^kLk{x) in Legcndrc orthogonal polynomials on the reference interval /, using the 
Cauchy-Schwarz inequality together with the properties |Lfe(±l)| = 1 and ||Lfc||g : = 2/(2k + 1). 
Another useful property will be a refined version of the classical inverse inequality ||w'||oj ^ 



P 2 

-77 IMIoj which holds for all v £ P p (I), and reads as follows: 



1/2 

v'\ki< (ly&K 1 I V«GP P (7). (2.6) 



A proof will be given in Sect. [6] 

Next, let d > 1 denote the space dimension of interest. Given d bounded intervals Jfc, 1 < k < d, of 
size iffc and carrying polynomials of degree pk, we form the corresponding d-tuples p = (pi, . . . ,pd) 
and H = (Hi, . . . , -ffd), we define the tensor-product space of polynomials of degree up to pk in the 
fc-th coordinate on the (hyper-) rectangle R — X fc=1 Ik as 

Q p (i?) = (g)P Pfc (4) . (2.7) 

fe=i 

Occasionally, we write p — p(R) and H = H(R). 

Given any integer < I < d, let J-'i(R) denote the set of all /-dimensional facets of R, i.e., the 
subsets ofdR obtained by freezing d—l coordinates to one of the endpoint values of the corresponding 
intervals. Thus, J^oiR) is the set of the 2 d vertices of R, Fd-i{R) is the set of all faces of R, whereas 
F d {R) = {R}. Obviously, if v € Q P (R) and F E Fi(R) for some 1 < I < d, then V\ F € Q p . (F), where 
p t = p*(F) is obtained from p by dropping the components corresponding to the frozen coordinates. 

Polynomials in Q P (R) are uniquely determined through their values at the tensorial LGL-grid in 

R 

d 

Q P {R) = X G Pk (h) = U = (£i,ji , 6, J2 , . . • , U,j d ) for < j k < p fc , 1 < k < d} , (2.8) 
fc=l 

where the £k,j k are the nodes of the LGL quadrature formula of order pk in Ik- Introducing the 
product weight = u>ij 1 W2j 2 ■ ■ ■ Wd,j d associated to each node £ G Q p (R), we obtain the quadrature 
formula 

^ v(Owt= I v(x)dx VveQ p (R). (2.9) 
Then, the discrete inner product in Q P (R) 



E u (0v(0n (2.io) 



defines a norm which, thanks to (2.4 1, is uniformly equivalent to the L -norm, i.e. 



IM|o,h < \Ho,r, p < (V3) d \\v\\ Q . R VveQ p (R). (2.11) 
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LGL-grids as well as discrete inner products and norms restrict to any facet of R in the obvious way, 



yielding objects with analogous properties (see Sect. 5.1 for more details) 



Finally applying the bound (2.5) and tensorizing in the remaining directions, we easily get the 



following "trace inequality", which will be crucial in the analysis of the DG discretization. 

Property 2.1. Let F € Fd-i{R) be a face of R, orthogonal to the k-th direction, for some 1 < k < d. 
Then, 

|| v || 0ii? <?*±l||„|| 0fl \fveq p (R). (2.12) 



2.2 Model problem and discretizations 

We formulate now various versions of what we call Spectral-Element Discontinuous- Galerkin (SE-DG) 
discretizations for second order boundary value problems. In principle, the developments presented 
below apply under fairly weak assumptions on the problem parameters. To keep the technical level 
of the exposition as low as possible we confine the discussion to the model problem 

-Au = f inn, u = on dn, (2.13) 

where C M. d , d > 2, is a bounded Lipschitz domain with piecewise smooth boundary and / 6 L 2 (il). 
Such domains can be partitioned into images of (hyper-) rectangles through smooth mappings (such 
as iso/sub-parametric mappings, or Gordon-Hall transforms). Again, in order to keep technicalities 
at a possibly low level, it suffices to treat unions of (hyper-)rectangles, i.e., we assume that the 
closure of £1 can be partitioned into the union of a finite collection TZ of (hyper-)rectangles R, which 
overlap at most through parts of their boundaries. Precisely, in this paper we will consider only 
geometrically conforming partitions TZ, characterized by the fact that any nonempty intersection 
between two elements R and R' in TZ is an Macet for both of them, for some < I < d — 1 . An 
example of a 2D conforming patch of elements is given in Fig. [I] It will be convenient to introduce 
the complex 

T x ■■= (J W) 

R.en 

of all Z-dimensional faces associated with the macro-mesh consisting of the elements R € TZ. 



Figure 1: A geometrically conforming patch of elements in M 2 , with the corresponding LGL-grids 



Recall from Sect. 



2.1 



that the fc-th component of the vector H = H(R) € Ri is the k-th side 
length of R. Likewise p = p(R) 6 N d denotes the vector of coordinatewise polynomial degrees of 
an element of Q p (R). The vector-valued, piecewise constant functions H and p are thus defined 
throughout f2 forming the function 8 — (H,p) which collects the approximation parameters for the 
chosen DG-SE discretization. 
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We assume the aspect ratios of the elements are uniformly bounded and that the distribution of 
polynomial degrees within each R as well as between adjacent elements R~ and R + is graded in the 
sense that 

max fe H k {R) ma,x k p k (R) p k (R ± ) / 01yl s 

; — - < 1, ; — - < 1, max max — ; < 1, R g 7c. (2.14) 

mm k H k {R) ~ min fc p fc (i?)~ ± k p k (RT) ~ v ; 

In particular, this implies that our meshes arc quasi-uniform. 

The approximation space for the DG-SE discretization is defined as 

Vs = {v g L 2 (f2) : VReTZ, there exists u fl g Q p (i?) such that = v R a.e.} ; (2.15) 

note that any v € Vs can be identified with the vector of polynomial functions (w_r)-Rg:r.- 

To introduce the standard notation for DG approximations of Problem (2.131, assume that F g 
J-d-i is an interior face shared by two elements which we denote by R~ and R + . If v is a piecewise 
polynomial function defined on both sides of F, which takes values v~ in R~ and v + in R + , then 
its jump [v]f across F and its average {v}f on i 7, are defined as 

[v] F = n fl - iir + n fl+)F uj^, {v} F = - (v |F + w+ ) , 

respectively, where ^ = — n^j+ ^ denote the unit normal vectors on F pointing to the exterior 
of R~ . Taking the homogeneous Dirichlet boundary condition into account, we set [v]p = hr,fV\f 
and {v}p = v\p when F C <9f2. 

Within this setting, we consider the Symmetric Interior- Penalty Discontinuous Galerkin Spectral- 
Element discretization of Problem (2.131 defined as (see OS]): find u g Vs such that 

a* s (u,v) = {f ,v) ,n V« g V 5 , (2.16) 
where the bilinear form on Vs x Vs is given by 

a* 5 (u, v)=^2 (Vu, Vv) 0> ij + X) ( ~ ({Vu} ' M)o > F ~ ({Vw} ' + TWf (M, H)o,f) ■ (2.17) 

The weights Wj? are defined as follows. Assume that the interface F is orthogonal to the fc-th 
direction. If F is shared by two elements R~ and R + , then we set 

f (p k (R-) + lf ( Pfc (fl+) + l) 2 \ 

where the particular choice of ujf will be justified later. On the other hand, if F is contained in dfl 
and belongs to the element R, we set 

_ {p k (R) + lf 

~ H k (R) ■ (2A9) 

Finally, the constant 7 > is chosen in order to guarantee the coercivity of the bilinear form as with 
respect to the DG-norm, defined in Vs by 

\M 2 dg,s= E l|V<* + 7 E ^IIHIIo,f- (2.20) 
Ren Fefd-i 

Indeed, the following proposition, which is a specialization of a well-known general result about DG 
methods to the present situation (see, e.g., |llj). confirms the uniform continuity and coercivity of 
the form a* Sl which implies well-posedness of the discretization scheme (2.16). 
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Proposition 2.2. There exists a constant 70 > such that for all 7 > 70 the bilinear form a| 

(2.21) 



defined in (2.171 satisfies 



a* s (v,v) ~\\v\\ 2 DG S VveVg. 
The constant 70 and the constants implied by the symbol ~ can be chosen independent of S. 
Proof. We just highlight the crucial step in the proof, which is the bound of each quantity 

|({V«}, [V]) , F \ < \\{VV}\\ ,F II M \\ ,F < ^ll{ Vu }Ho,F + ^|| M f ,F, 

where r\ > is a suitable weight allowing us to control these terms by the two terms appearing in 
the DG-norm. Assume first that F is an interior face, orthogonal to the k-th direction and shared 
by the two elements R~ and R + . Then, denoting by the partial derivative with respect to the 
k-th. coordinate, we obtain 

II{V4IIo,f = \\\(d k v-) ]F + (d k v+) lF \\l F < ]T mv^Wlp < ^^119^11^ , 



where we have used the trace inequality (2.121. A similar estimate holds when F is a boundary face. 
The rest of the proof is standard. □ 



2.3 Numerical integration 

Computational efficiency is usually enhanced, within a nodal-basis approach such as the one adopted 
here, by inserting high-order numerical integration (NI) into the Galerkin scheme at hand. This 
is accomplished by replacing the exact L 2 inner products on any element R and any face F by 
suitable discrete inner products based on Legendre-Gaufi-Lobatto quadrature formulas. While on 



each element R it is enough to use the inner product defined in (2.101, some care is needed on a face 
F, where restrictions of polynomials of different degrees may come into play. Let us detail the latter 
situation. 

Assume that F € J^d-i is shared by two contiguous elements R , and is orthogonal to the k-th 
direction. Let pf- be the vectors of (d— 1) polynomial degrees, obtained from p ± = p(R ± ) by deleting 
the k-th entries. If v € Vs, and if iF denotes its restriction to R ± , then (v ± )^p belongs to Q p ±(F). 
Let us introduce the vector = max(pj,p+), where the maximum is taken component- wise. Then 
both (i> - )|f and (f + )|F belong to Q Pr (F). Thus, it makes sense to consider the (d— l)-dimensional 
LGL-grid Q Vr (F) on F corresponding to the polynomial-degree vector p* , as well as the corresponding 
quadrature nodes. The induced inner product in Q Pi (F) will be denoted by 

= u(C)v(C)w c . (2.22) 

CeS P „(F) 

With these notation at hand, we assume that / 6 C°(f2) and we obtain the modified DG-NI (Dis- 
continuous Galerkin with Numerical Integration) scheme: find u G Vg such that 

aS,Ni(«.«) = (/,«)n,Ni Vv€V s , (2.23) 

where 

^,Nl( M ' v ) = ^2 ( Vu > ^ v )o,R,p 

Ren 

^2 ( ~ ({ Vu }' M)o,F,P» ~ [ u ])o,F,p, + 7^f(M, M)o,F,p.) , 



(2.24) 
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and 

E(/>v)o,h*. (2-25) 
Ren 

Note that at the algebraic level, the computation of the interface terms may require suitable transfer 
operators between LGL-grids. 

Proposition 2.2 admits a completely analogous statement for the form ag NI , which can be estab- 



lished through repeated applications of the norm equivalence (2.11 1 in each R and F. 
2.4 The SE-DG reduced bilinear form 



In view of Proposition 2.2 and its Nl-counterpart, both forms a* s and a^ NI are uniformly equivalent 
to one of the reduced bilinear forms 

as(u, v)=^T (Vu, Vv) 0tR + 7 

u F (HM)o,f, (2.26) 

Ren Fe^d-i 

or 

as,m(u, v) = (Vw, Vw) ,_r, p + 7 ^ ^(M) [ v ])o,f,p, ■ (2.27) 
Ren FeJi-i 

Therefore, we will continue our discussion by investigating efficient preconditioners for these forms. 



In particular, theoretical considerations will be developed for the conceptually simpler form (2.261 



although numerical results will be obtained using the complete, yet computationally more efficient 



form (2.241 



3 The auxiliary space method 



A suitable conceptual framework for developing preconditioners for the linear system (2.16) is offered 
by the so called auxiliary space method (see, e.g., [3 |2SJ [T2])- For convenience of the reader we 
briefly recall now the general setting adopting the abstract framework, see [2UJ HI] . 

Let V be a finite-dimensional space contained in some functional space X, and let a : V x V — > K 
be a given symmetric positive-definite bilinear form. Let us introduce an auxiliary space V, still 
contained in X and finite-dimensional, which carries a symmetric positive-definite bilinear form 
a : V x V — > M. Let us assume that two symmetric positive-definite bilinear forms are defined on 
the sum V = V + V, say a, b : V x V — > M, which satisfy the following conditions: 

i) a is a spectrally equivalent extension of both a and a, i.e., 

a(v, v) ~ &(v, v) Vw G V , (3-1) 

a(v, v) ~ a(v, v) e V . (3.2) 

ii) b dominates a on V, i.e., 

a(v,v) <b(v,v) VveV; (3.3) 

iii) there exist linear operators Q : V — > V and Q : V — > V such that 

a(Qv, Qv) < a(v, v) e V , (3.4) 

a(Qv, Qv) < a(v, v) £ V , (3.5) 

and 

b(v-Qv,v-Qv) <a(v,v) VveV, (3.6) 

b(v-Qv,v~Qv) <a(v,v) Vu e V . (3.7) 
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Remark 3.1. If condition ii) above is replaced with the stronger condition: 
ii') b dominates a on V, i.e., 



a(v,v) < b(v,v) VveV, (3.8) 



then one can skip checking inequalities (3.4 1 and (3.5 1 in Hi). 



Proposition 3.2. Under the assumptions (3.1 1— (3.7), the following norm equivalence in V holds: 
ca(v,v) < inf _ {b(w, w) + a(v, v)} < Ca(v,v) VveV, (3.9) 

w e V, v e V 

v = w + Qv 

where the constants c and C only depend on the constants implied by the assumptions (see for 
their expression) . 

Remark 3.3. It is worth mentioning the special case when V can be chosen as a proper subspace 
of V . Then V — V and natural choices for a and a are the form a itself, whereas an obvious choice 
for Q is the canonical injection. In such a situation, the only assumptions that need be checked for 
Proposition \3.2\ to hold are: 

a(v,v) < b(v,v) VveV (3.10) 
and the existence of a linear operator Q : V —> V such that 

b(v - Qv, v - Qv) < a(v, v) Vw e V . (3.11) 



At the algebraic level, Proposition 3.2 has the following main consequence. Let A, A and B 
denote the Gramian matrices for the bilinear forms a, a and b (restricted to V x V) with respect 
to suitable bases of the spaces V and V. Let S be the matrix representing the action of Q in these 
bases. 



Corollary 3.4. Under the assumptions ( 3.1 ) — ( 3.7 1 , let Cb and be symmetric preconditioners 
for B and A, respectively, satisfying the following spectral bounds: 

Amax(CBB), A max (C^A) < A max , A m j n (CBB), A m i n (C^A) > A m ; n . 

Then, Ca = Cb + SC^S T is a symmetric preconditioner for A, and the spectral condition number 
of CaA satisfies 



Miiax 



k(CaA) < 



4 Orientation 

We shall employ the auxiliary space method for preconditioning the linear system arising from the 
SE-DG - Spectral- Element Discontinuous-Galerkin discretization introduced in Sect.[2j 

Perhaps a few comments on the possible choices of auxiliary spaces are in order. It has been seen 
in |11) 112] that for bounded polynomial degrees a low order conforming finite element subspace works 
well. Rougly speaking, it then sufficies to apply simple relaxation to the "nonconforming" part of 
the DG finite elements space and to employ multilevel preconditioning techniques to the conforming 
low order subspace. Since constants arising in this analysis depend on the polynomial degree in 
the DG-space this strategy cannot be expected to work in its original form in the present context. 
Instead one learns from experiences in the spectral element method (see e.g. [T5J HH H3] ) that 
suitable auxiliary spaces should be based on low order finite elements on LGL-grids, introduced in 
Sect. EQ 
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However, then several serious difficulties arise. The first one is rather obvious, namely how to 
precondition in the case of very large degrees corresponding low order finite element discretizations 
since the particular associated grids are not nested and typical multilevel techniques cannot be 
applied directly? Moreover, multilevel techniques work well, as mentioned above, when dealing with 
(globally) conforming discretizations. This latter issue is more subtle but turns out to pose serious 
obstructions as soon as one uses different polynomial degrees on different patches. Specifically, LGL- 
grids on adjacent patches with different polynomial degrees do not match and globally conforming 
subspaces are too thin for preconditioning. Therefore, it seems that an additional "second stage" 
auxiliary subspace is needed. For this purpose we shall construct suitable low order conforming 
finite element subspaces on certain dyadically refined (but anisotropic) grids associated with the 
LGL-grids, see Sect. [6] giving rise to the DFE-CG - Dyadic-Finite-Element Continuous- Galerkin 
schemes. These grids can be nested and thus used to form sufficiently "rich" globally conforming 
subspaces to which multilevel preconditioners apply. 

Now the question arises "why not using the latter spaces directly as auxiliary spaces for the SE- 
DG discretizations?" Again, a second look at this option reveals that the jump terms in the DG 
bilinear form cause serious problems when trying to fulfill the above requirements in the auxiliary 
space framework. In fact, although the degrees of two adjacent polynomial patches may differ the 
jump across the interface could vanish since one of the traces may happen to degenerate to the 
common lower degree. However, since the corresponding LGL-grids do not match the respective 
associated dyadic grids differ as well and will therefore give rise in general to nontrivial jumps of 
the corresponding piece wise linears. This incompatibility of the jump terms seems to hinder the 
verification of the ASM conditions. 

Therefore, it seems that one has to resort to an iterated application of the auxiliary space method. 
The preceeding remarks concerning the compatibility of jump terms in different discretizations sug- 
gest to reduce the problem as quickly as possible to a conforming one for which one can employ 
standard energy inner products. Therefore, as a first-stage auxiliary space we choose the largest 
conforming subspace of V$, referred to as SE-CG - Spectral- Element Continuous- Galerkin space. 
This is already interesting in its own right since it can be combined with (nearly optimal) domain 
decomposition preconditioners, see |16j . 

Since local FE spaces on contiguous elements with LGL-grids can be glued together to produce 
non-trivial globally conforming spaces only when the polynomial degrees are all equal, we shall use 
this concept only implicitly and propose to pass from SE-CG directly to DFE-CG. Once we have 
arrived at that stage, one can proceed to preconditioning this latter problem by suitable multilevel 
techniques. Due to the anisotropies of the involved dyadic grids this is, however, not quite straight 
forward. We develop therefore a multi-wavelet based preconditioner whose detailed elaboration is, 
however, deferred to the forthcoming part II of this work. 

In conclusion, we propose the following route to an iterated auxiliary space preconditioner: 

• IAS-PATH: SE-DG -4 SE-CG -4 DFE-CG . 

The first stage will be discussed in Sect. [5] whereas the second one will be detailed in Sect. [7j 

We proceed collecting the necessary technical prerequisites for the LGL-grids in Sect. |5.1| and for 
the associated dyadic grids in Sect. [6] 

5 Stage SE-DG -+ SE-CG 

In this section, we investigate the use of conforming spectral elements as auxiliary spaces. Before 
defining the constitutive ingredients and checking the ASM assumptions, we introduce some addi- 
tional notation and we add a few remarks, which will be helpful here and throughout the remainder 
of the paper. 
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5.1 Auxiliary material 

Given two integers < I, m < d and a facet F G Fi, we set 

F m {F) = {GeF m : G C F if m < I , G = F if m = I , G D F if m > 1} . (5.1) 

In particular, J-'i(R) has been already denned as the set of all Z-facets of an element R G 1Z, whereas 
1Z(F) = Fd{F) is the set of all elements R containing the facet F. 

Given any R G 1Z and the corresponding LGL-grid G P (R), if F G Ti{R) for some 1 < I < d — 1, 
then Gp(R) O F is the LGL-grid of the induced quadrature formula on F, which will be denoted 
by Q p (F,R). The order of the formula, denoted by = p(F,R) G N l , is obtained from the vector 
p = p(R) by deleting the d— I components corresponding to the frozen coordinates of F. The weights 

F R 

of the formula, denoted by We ' , are connected to the weights of the original formula on R by 
the relation 

W( = Wfi ■■■w fd _ l wf R V££Gp(F : R), (5.2) 

where each uif i is a boundary weight of the univariate LGL formula along one of the frozen coordi- 
nates of F. This implies that if the facet F' belongs to Fi-i(F) C J r /_i(i?), then 

wf > R = w f wf' R V£ G G P (F', R) , (5.3) 

for some boundary weight Wf. 

Next, let us draw some consequences from the assumptions (2.14) of quasi-uniformity of the mesh 
and polynomial degree grading. Under these assumptions, it is easily seen that given any R G 1Z, 
there exists a real H = H/j > and an integer p = pjj > 0, such that for all R! 6 1Z satisfying 
R' n R ^ 0, it holds 

H k (R')~H R , p k (R')~p R , l<k<d, (5.4) 

i.e., lengths and degrees are locally comparable. More generally, for any face F E Fi, I < d, one 
can define analogous "representative parameters" Hf,Pf, for instance, by averaging corresponding 
parameters from intersecting faces. This implies, on the one hand, that each boundary weight wj, 
j € {0,pk(R')} of any univariate LGL formula used in the definition of the tensorial LGL formula 
on R' satisfies 

H , , 

Wj ~ , (5.5) 

and, on the other hand, that each face F <E Fd-i(R') carries a DG weight cjp satisfying 

P 2 

uj f ~ — . (5.6) 

In particular, all the faces F G Fd-i having a nonempty intersection with an element R carry weights 
cop of comparable size. 



5.2 Definition of the ASM ingredients 

Referring to the notation of Sect. [3j let us choose as V the space V$ introduced in (2.151, whereas 

applies, and we are only 
V — > V for which 



3.3 



as V we choose the space Vg = Vg PI Hq(Q). Since V C V, Remark 
required to define two bilinear forms a, b : V x V R and a linear operator Q 
conditions (3.101 and (3.111 are fulfilled. 



The bilinear form a will be the reduced form a§ defined in (2.26). 
the standard inner product in Hq(CI), since 



Note that on Vg it reduces to 



as(u, v) = 2J (Vu, Vw) 0: fl = (Vw, Vw) 0: o Vw, v G V s 
Ren 
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The definition of the form b is inspired by the inverse inequality (2.6) and requires the following 
notation. Let R E 7Z be any element, and £ E G P (R) be any LGL quadrature node in R, with 
associated weight = Wj 1 Wj 2 ■■■Wj d ; we denote by w^ t k the factor Wj k coming from the k-ih 
direction. Then, we introduce the weights 

w * = ( e w ^ ) w * = e w ii ( n w « d) • ( 5 - ? ) 

\fc=l / k=l j^k 

We also introduce strictly positive coefficients c^, satisfying the requirement ~ 1 (precisely, they 
should be bounded from above and from below independently of £, p and H under the assumptions 
(2.14)). They will be chosen with the aim of enhancing the effects of the ASM preconditioner; we 
refer to Sect. [8] for the details. 

Definition 5.1. The bilinear form b$ : V$ X Vs — > K is defined as follows: 

bs(u,v) = ^2 M w ' v ) > with b R (u,v) := u(f) v(£) c e W € . 

Note that the bilinear form 65 is defined strictly elementwise, i.e. it does not involve any coupling 
between different elements R, which is essential for the efficient applicability of the preconditioner 
B in Corollary |Q| 

Before giving the definition of the operator Q$ : Vs — > Vg, we introduce a new assumption on the 
local distribution of the polynomial degrees, which indeed poses a restriction only for d > 3. 

Assumption 5.2. For 2 < I < d — 1 and for each F 6 Ti, there exists R E TZ(F) such that 

p(F, R) < p{F, R') VR' E K(F) 
(the inequality being taken componentwise). 

Note that the property is trivially true for I = 1 since p(F, R) E N, and for I = d since 1Z(F) = {F}. 
Thus, for each F E J~i with 1 < I < d, we are entitled to select, once and for all, an clement R^(F) 
such that 

R*(F) := argmin ra(i , )P (F, It), ffl(F) := p{F, R*(F)). (5.8) 

The LGL-grid g p (F,R i (F)) will be denoted by G p t(F). Finally, for each vertex x E F of the 
decomposition we select an element R^(x). 

The definition of the operator Qs will be accomplished through a recursive procedure, as follows. 

Definition 5.3. For < I < d, set 

Fi = U F ■ 

Given any v E Vs, let us define the sequence of piecewise polynomial functions qi : Fi — >• R by the 
following recursion: 

i) For any x E Fq, set qo(x) = if x E d£l and qo(x) = v^t^^x) if x E f2. 

ii) For I = 1, . . . , d, define qi on Fi as follows: for any F E T\, set qi\p = if F C dfl, otherwise 
define qi\p by the conditions 

Qi\f G Qpt(F) such that 

«|f(0="*»(f)(0 ^^G P »(F)\diF , (5.9) 
qi\ F (0 - qi-i(0 Ve € g p t(F) n d t F , 

where diF denotes the boundary of F as an l-dimensional manifold. 



Multilevel Preconditioning of DG-SEM, Part I: Geometrically Conforming Meshes 



13 



Finally, we set 

Qsv = q_d ■ 

Remark 5.4. The above recursion defines a linear operator Qs from Vs into Vs- 



(5.10) 



In fact, the linearity of Qs is obvious. Moreover, note that Qsv £ Vs since p*(R) = p(R) for all 
R £ 1Z. Furthermore, one has the property 



%Fi- x - qi-i 



1 < l< d . 



(5.11) 



Indeed, if F G T\ andF' G Jj-i(-F), theng ; | F / G Q p (f>,rs(f)){F') whereas qi-x\F> G Q p (f> W(f> ))( f ')- 
The minimality property (5.8) applied to F', together with the last set of conditions in (5.9 1, imply 
that qi\pi — qi-i\F'i whence (5.11 ) follows. This property implies that Qsv is continuous throughout 



f2, and vanishes on dQ.. We conclude that Qsv belongs to Vs- 



5.3 Check of the ASM assumptions 



We shall now first verify (3.101 in Remark 3.3 To this end, the following immediate consequences 



of (5.51 and (5.6 1 concerning the weights introduced in (5.7 1 will be useful. 

Lemma 5.5. Let R £lZ be arbitrary. For any face F £ Td-\ o,nd any £ £ Q P (F,R), one has 

~ u>fw^' R ■ 

Proposition 5.6. For bs given in Definition \5. 1\ one has 

as(v,v) < bs(v,v) Vu G Vs . 



Proof. For any R G 7Z, using (2.3 1, ( |2.4[ ) and (2.6) together with tensorization, one obtains 

\vv\\i R = ih\\d k v\\i R <j2 E <o 2w il(U w ^)= E «(0 2 we. 



fc=i 



whence 



E iiv«ii§,« 



On the other hand, let F £ Td-\ be a face shared by two elements, say R* 1 . Then, recalling the 
definition of pf given in Sect. |2.3l we have 



ii m \\1,f < E ii^iis.* < E n^iio,^ = E E V fl± 

± ± ± ?6e P (-F,fl±) 

Multiplying by w F and using Lemma |5. 5 1 for both R = R ± , we obtain 



^iimiIo,f<E E (« ± (0) 2 ^ 5 ± <E E ( v± (0) 2 w± 

± cee P (F,fl±) ± 5ee P (-R±) 

A similar result holds for the faces F £ Td-\ sitting on the boundary of f2. Hence, 

E wf 'I M Ho,-f ~ 6 «( u > u ) » 



and the result is proven. 



□ 
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We now focus on verifying assumption (3.111. This will be accomplished through a sequence of 
intermediate results. 

Lemma 5.7. For any v € Vg, let its restriction to R S 1Z be denoted by vr. The following bound 
holds 

bs(v - Q s v,v - Q s v) < J~] uj F \\v R - q d -i\\l , F , v &V s , 

Fefd-i Ren(F) 



where qd-i is constructed in Definition 5.3 in order to define Qsv. 
Proof. By definition of b$, one has 

b s (v-Q s v,v-Q s v)=J2 E \( v R-Qsy)(0\ 2 c 6 W^ 
Ren £eg P (R) 



In view of (5.10 ), (Q$v)\r, — qd\R and q_d coincides with vr on Q p (R)\dR and with qd-i on Q p (R)f]dR. 
Thus, since ~ 1 

b s (v - Q s v,v - Q s v) ~ £ E l(«H-5«i-i)(0| 2 Wj 

Ren £eg P (R.)ndR 

^ E E E i(«*-&-i)(0i a w« 
^ E ^ E E k^-^-x)^! 2 ^, 

F6fd-i Ren(F) £eg p (F,R.) 



where we have used Lemma 



5.5 



in the last inequality. Finally, we observe that vr\f € Q. p (f,r)(F), 
whereas q d _ 1{F € Q ptt (F) (F) <T0p(F„R)(F). Thus, ( |2.11[ ) yields 

E K^-^-iXOI 2 "^- ll u -R _ ^-illo,F 

proving the assertion. □ 



Lemma 5.8. Let v £ Vs be arbitrary, and let qi, < / < d, be the sequence built in Definition \5.3\ 
in order to define Qsv. For 1 < I < d— 1, given any F £ T\ and any R € 1Z(F), the following bound 
holds 

\\vR-qi\\l,F< E II Mg 11^ + 4 E IK«(F)-«-illo,F', (5.12) 

Gefd-^F) P f'EJi-ifF) 

where [v]g denotes the jump of v across the internal face G, or v\a when G C dfl. Here, H = 
and p = pf are representative values of the local meshsizes and polynomial degrees, introduced in 
(5.4 1. On the other hand, for any vertex x S J-q and any R € TZ(x), one has 

\(v R -qo)(x)\ 2 < E IHg(^)| 2 - (5-13) 

Proof. Assume first I > 0. If F C dR n dtt, then 

\\vR-qi\\l F =\\v R \\l F < E HMgIIo,f> 

Ger d -i(F) 

which is a particular instance of ( |5.12| . Let us then assume F <£_ d£l and, for simplicity, let us set 
i?" = R^(F). Observe that the definition of qi\ F on Q p t{F) n d[F given in ( 5 .9 1 can be rephrased as 
Qi\f(0 = v R f(£) + (qi-i(0 - vr»(0), or equivalently, 

qi\F = v R t\ F + E (Qi-i{Q- v Rt(Q)i>( , 

?ee »(F)n9,F 
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where ip^ G Q p t (F) denotes the Lagrangean function associated with the node £ of the LGL-grid on 
F of order pK Thus, using once more ( 2.11| , we obtain 

\\vR-qi\\l F <\\vB -v Ri \\l F + £ - VR«(0\ 2 ™f' R> =-A 2 + B 2 . 

As for the first quantity, there exists a sequence R l G 71(F), i = 0, . . . ,m < d — I such that R° = R, 
R rn = i?H and i?" 1 flff e F d -i(F) for i = 1, . . . , m. Thus, 

a 2 <EIK*-^-*IIo,f< J2 WMoWIf- 

*=1 GeJ"d-x{F) 



On the other hand, using ( |5.3[ ) and (5.5 1, we have 
B 2 < E E \Ql-i(0-vm(0M' R * 



< 



H 



E E kvi(e)-^«(0i 2 ^. 



2 

s 



E H w -R» -«-illo,F' . 



F'e^i-i(-F) CGS P (F',_R») 



where the last equivalence follows, as in the proof of the previous lemma, by observing that Vp }\ F ' G 
%(f<,r*){F') whereas q t _ llF , G Q p t {Fn (F') C Q p(F%W) ( J F 1 ')- Thus, (J532J is proven. Finally, fl5T3t 
is trivial. □ 

Lemma 5.9. Lei F e Fi for some < I < d — 2, and /e£ G € Fd-i(F). Define p* — max ^p(R) 
(componentwise). Then, 



d-l-l 



\Hl F <\\v\\l G WeQ p .(G) 



Proof. Using several times (|5.3| and (|5.5|), as well as (|2.1l|, we have 

,.G 



\l G * e i««)i a «f * E noM 

€e6p-(G) ?eS p .(G)nF 

" ' 1 E woiV 

?eS P *( F ) 



H 



d-l-l 



\Q,F 



□ 



We are now ready to prove condition (3.11 1. 



Proposition 5.10. For bg given in Definition 5.1 and for Q$ given in Definition 5.3 one has 

bs(v - Q&v,v - Q$v) < a s (v,v) Vw G V s . 



Proof. First, we observe that the cardinality of any set F m (F) defined in (5.1l is bounded by a 
quantity depending only on the dimension d. We start from the bound given by Lemma |5.7| and 
focus on any face F G Fd-i- Using inequality (5.12) of Lemma [5. 8| with I = d — 1, we get 

E \\ V R - 9d-l\\o,F £ II M llo,F+ ~2 E W V RHF) -Qd-2\\l,F> ■ 
R.e1Z(F) P F'€J r d _ 2 (F) 
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A further application of Lemma |5.8| to each summand on the right hand side of the above inequality, 
now with I = d — 2, yields 



E \\ V R ~ Qd-l\\l,F ^ 
Reiz(F) 



||2 , _H 

I 0,F ' n 
P 



E E 



Ma Wo.f' 



E E W V RHF') -Qd-sWl,! 

F>£F d ~ 2 (F) F"£F d - 3 (F') 



5.9 



with I — d — 2 yields -^|| [w]g ||q f> ~ II H Ho G- At this P om t, let us introduce the set 
n F ^ 0} of all faces intersecting the face F, and 

E E iiMg||oV< E hmiIo.g 



Lemma 

FjL^F) = {G 6 J^d-i '■ G n F ^ 0} of all faces intersecting the face F, and let us observe that 
H 



F'<£F d - 2 (F) G<£F d -!(F) 

On the other hand, we have 



E E i 

F'eF d „ 2 (F) F"£F d _ 3 (Fi) 



VRi(F') - Qd-3\\o,F" ~ 



Thus, 



Y IK-%-lHo,F< E HHIIo.G 
ReTZ(F) GeF^iF) 



GeF^F) 



E E W V R~ Qd-3\\l,F" 

F"eF d - 3 {F) R&K{F") 



E E ll U « -^-3|lo,F" 

F"eF<i_3(F) Ren(F") 



We now proceed recursively, using Lemmas 5.8 and 5.9 with I = d — 3, d — 4, At the j-th stage 

of recursion, we obtain 



E ll w R-<ta-i|lo,F ^ 



E II H 



2 

lo.G 



H 



3-1 



R£1Z(F) 

When j — d, we use ( 5.13| and ( |5.5| to finally get 

E W V R- Vd-l\\l :F ^ E 



E E ik-^-jIIiIf' 



Mllo.G 



Ren(F) 



At last, we observe that ~ ujr for all G € J r £ J_ 1 (.F) by (5.6 1, so that, going back to Lemma 
we conclude that 

bs(v - Qgv,v - Q s v) < Y w g||HIIo,g < a s (v,v) . 



5.7 



□ 



5.4 Summary of Stage I 

Thus, all requirements of the Auxiliary Space Method are satisfied for the quantities defined in 



Sect. 5.2 Note that the Gramian B of the form bg is diagonal so that Cb = B~ : is efficiently 
applicable and trivially satisfies the respective requirements in Corollary |3.4| The quantitative per- 
formance of this preconditioning stage will be discussed in Sect. [8] So it remains to construct a 
suitable preconditioner for the "auxiliary space problem" involving the bilinear form a. Several 
possible routes to proceed from here suggest themselves. The first one is employ a domain decompo- 
sition strategy in the spirit of [25j . This is carried out in [16] where a dual-primal preconditioner for 
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the auxiliary problem is formulated and analyzed. Specifically, the resulting condition numbers are 
shown to exhibit only a mild square logarithmic growth with respect to polynomial degrees similar 
to those obtained in [25] for constant degrees. 

Here we pursue a different route trying to develop a preconditioner that eventually gives rise 
to uniformly bounded condition numbers for the auxiliary problem and hence also for the original 
problem. Recall that the main obstruction to employing the low order LGL-grid finite elements as 
a next-stage auxiliary space for the conforming high order space is the non-matching of the LGL 
grids of different degrees at element interfaces which, in turn, is closely related to the fact that the 
LGL-grids are not nested. A key element of our approach is to construct certain low order auxiliary 
spaces that are in some sense "stable companions" of the low order LGL finite elements but that are 
nested. However, having overcome the principal hurdle of non-nestedness certain further difficulties 
arise that will be addressed throughout the remainder of this paper. 

6 Piecewise polynomial spaces over LGL and dyadic meshes 

This section contains preparatory material for the realization of the subsequent stage SE-CG — > 
DFE-CG, see Sect. [4] After recalling some further properties of the LGL-grids introduced in 
Sect. [2] we turn to the following crucial task. We develop an algorithm to associate to any LGL- 
grid a dyadic grid in such a way that nestedness is guaranteed as the polynomial degree increases. 
Several approximation operators related to these grids are defined, and certain stability properties are 
established. Finally, the previous univariate results are extended by tensorization to the multivariate 
setting which eventually allow us to use these dyadic grids in place of LGL-grids for the construction 
of low order auxiliary spaces. 

Let us introduce a few notations which will be used throughout the section. Given a bounded 
interval / = [a, b] C R, consider an ordered grid Q — {£j : < j < p] C /, with 

a = Co < ■ ■ ■ < < & < • • • < £p = b . 

The collection of intervals Ij = 1 < j < p, defines a finite partition (or mesh) of /, which 

will be indicated by T = T{Q). Conversely, any finite partition T of I comes from a unique ordered 
grid Q = G(T); for this reason, the two concepts will be often used equivalcntly in the sequel. We 
denote by hj = £j — = the length of the interval Ij 6 T. Furthermore, we associate to the 
partition T the space of continuous, piecewise-linear functions 

V h (T) = {ve C°(I) : v\ I} € Pi , VIj € T} . (6.1) 

The following definitions will be crucial in the sequel. 

Definition 6.1. An ordered grid Q = {£j : < j < p} in I is locally C g -quasiuniform if there exists 
a constant C g > 1 such that 

C; l <^-<C g , l<j<p-l. (6.2) 

Definition 6.2. An ordered grid Q in I is locally (A, B) -uniformly equivalent to another ordered 
grid Q if there exist constants < A < B such that 

VI 3 eT(G), VIieT(G), IjDI^V) =► A < Hq- < B . (6.3) 

W 

Examples are given by the LGL-grids and the associated dyadic grids, as discussed below. 
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6.1 LGL-grids and finite elements 

Given any integer p > 0, let G P (I) be the LGL-grid of order p in / and let T P (I) — T(Q P (I)) be 
the corresponding partition. We recall that the points in G P (I) are symmetrically placed around the 
center of the interval, and the intervals in T P (I) are monotone in the sense that 

I 3 ,I ]+1 eT p (I), Ij, C [a,(a + b)/2] => \I J \<\I J+1 \, (6.4) 

with the analogous reverse relation for the right half [(a + b)/2, b) of I. Furthermore, the variation 
of the interval lengths is uniformly smooth, in the sense that G P (I) is locally C g -quasiuniform for a 
constant C g > 1 independent of p. These results can be found, e.g., in |5]. 

Another relevant property concerns the locally uniform equivalence of LGL-grids of comparable 
order; we refer to [S] for the proof. 

Property 6.3. Assume that cp < q < p for some fixed constant c > 0. Then, Qq(I) and G P (I) 
are locally (A, B)- uniformly equivalent, with A and B depending on the proportionality factor c but 
independent of q and p. 



A crucial role in what follows will be played by the relation between the high-order polynomial 
space V P (I) := P p (/) and its low-order companion finite element space Vh, p (I) = Vh(Tp(I)) of the 
continuous, piecewise-linear functions on the partition Tp(I). Let us introduce the interpolation 
operators at the nodes of G P (I) 



l p :C°(I)^V p (I) , 

fot;)&)=«&), 



l htP : C°(I) -+ V htP {I) , 



0<j< P , Vve C°(I) 



(6.5) 



Then, v i— > U/, = Xh, P v defines an isomorphism between P p (J) and V/j lP (J), with inverse »), 4 u = 
I p Vh- The next property provides uniform equivalences between norms of global polynomials and 
their piecewise affine interpolants at the LGL-grid, and vice-versa, which are at the core of the 
subsequent applications to preconditioning (see [131 [23] ) • 



Property 6.4. One has 



and 



|M| 0i/ ~ K||o,j WeP p (/) 



IKIk/^KIIo,/ v ue p p (J) 



(6.6) 
(6.7) 



where the involved constants are independent of p and H = b — a. 



Note that taking as v in (6.6) each Lagrange basis function at the LGL nodes and using (2.2 1 
and (6.2 1, it is easily seen that the size of each interval Ij is comparable to that of the LGL weight 
associated with any of its endpoints, in the sense that the following bounds hold, uniformly in p and 
H: 



1 mm 

i<i<p Wj 



< max — <■ 1 . 



l<j<p w 



(6.8) 



As a consequence, (6.7 1 together with (6.8) and (2.2) provide a simple proof of the inverse inequality 
(2.6 1. Indeed, one has 

IKIlo,/ 



E 

3=1 



V(Co)+E 

3 = 1 
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6.2 Dyadic meshes 

The subsequent stage SE-CG — > DFE-CG will rely upon the use of auxiliary spaces based on 
dyadic grids that are in a certain sense associated with LGL-grids, but have the additional advantage 
of giving rise to nested conforming spaces. Their construction, based on recursion algorithms, is 
detailed in [5] (see also [S] for an extensive quantitative investigation of their properties). Hereafter, 
we summarize the main results that are needed in the subsequent analysis, referring to [S] and [5] 
for their proofs. 

We wish to construct for a given ordered grid Q in / an "equivalent" (in the sense of the definition 
above) dyadic mesh T>. To this end, it will be convenient to introduce for any interval D 

1(D,G) :=argmax{|J,-| : I s £ T(Q), Dn/^H} . (6.9) 

Given a real a > and an initial dyadic partition Uq, a dyadic partition T> can be constructed 
recursively as follows: 



Dyadic [g,V ,a] V: 
(i) Set V:=V Q . 

(ii) While there exists D € T> such that 

\D\>a\I{D,g)\ (6.10) 
split D by halving it and replace it by its two children D' , D" , D = D' U D" , i.e. 

(V\{D})U{D',D"}^V. 



Such a dyadic mesh generator, applied to the sequence of LGL-grids g p — g p (I), produces meshes 
T>* := Dyadic [g p , {I}, a] which enjoy additional structural properties useful for our purposes. First, 
as the sets g p are symmetric around the center of the interval, so will be the grids T>*, always con- 
taining the center as a node. Next, a major requirement on the dyadic partitions associated with the 
g p 's is that they are nested, which will later be essential in the context of the ASM. Unfortunately, 
numerical evidence shows that the grids T>* are not necessarily nested, although exceptions seem to 
be very rare (see [S]). Thus, in order to ensure nestedness of the generated dyadic meshes, we shall 
employ the following recursive definition: 



y p " Ly P ' 

(i) Given a > 0, set 
(ii) For p > 1, given X> p _i, set 



V 1 : =Dyadic[0i, {/},«]. 



U p := Dyadic[^ p ,X> p _i,a]. (6.11) 



Property 6.5. For all p > 1, the dyadic meshes T> p are nested and locally (A, B) -uniformly equiv- 
alent to g p with constants A, B specified as follows: 

I 7". | Of 

VDen pi yi,eT P , I, HD ^0 => A ■= or 1 < j-^j < 9 _ 1 - =: B . (6.12) 

Furthermore, 

cardPp ~ cardCJp . (6.13) 
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In order to keep the cardinality of V p as close as possible to that of Q p , numerical evidence suggests 
to choose the parameter a > 1, to guarantee gradedness of the dyadic mesh, but close to 1, to avoid 
excessively many points in the dyadic mesh, see [3] for more details. 

The following useful generalization is an easy consequence of Properties |6.3| and |6.5| 

Corollary 6.6. Assume that cp < q < p for some fixed constant c > 0. Then, T> q is locally (A, B)- 
uniformly equivalent to both Q p and T> p , with A and B depending on the proportionality factor c but 
not on p. 



6.3 Stability results for univariate interpolation operators 

In the rest of the paper, we assume that we have fixed a value a > 1 close to 1 and consequently 
for any interval / we have generated by the recursion (6.111 the sequence T> p (I) of dyadic meshes 
associated with the LGL-grids G P (I)- Let Gd,p(I) denote the set of dyadic nodes in I which define 
the partition T> p (I), i.e., such that T) p (L) = T(Gd,p(I))- Furthermore, let us introduce the finite 
element space 

V h , D , p (I) = V h {V p {I)) = {ve C°(I) : v\ D ePi,VDe V p {I)} (6.14) 
of the continuous, piece wise linear functions on the partition T> p (I), and let 

lh,D, P ■ C°(I) -> V hjDjP (I) , (l hjD , p v)(C) = v(() VCe£ Ap (J), VveC°(I), (6.15) 

be the interpolation operator at the nodes of Gd. p {I)- We also introduce the composite interpolation 
operator 

K>h,D, P ■ C°(I) V h ,D, P (I) , K-h,D, P = Zh,D, P °Zh, P ■ (6.16) 
We proceed now to establish various stability results involving these operators. 

Lemma 6.7. Let Q be any ordered grid in I, which defines a partition T = T{G), and let Xg : 
H 1 (I) — > Vh(T) be the associated piecewise linear interpolation operator. Then, 

MM<Hy VveH\i). 

Proof. Let Q = {£,• : < j < p}, with hj = £j Then, 

j—l^ij-l V 3 / j = i 3 

Writing v(£j) — «(£j-i) = J^ 3 1 ^ and using the Cauchy-Schwarz inequality, we obtain 



C-i ^ dx 



whence the result follows. □ 



Lemma 6.8. Let Q and Q be ordered grids in I , with associated partitions T = T(G) and T = T(G)- 
Assume that Q and Q are locally (A, B)-uniformly equivalent. LfTg : H 1 ^) — > Vh{T) is the piecewise 
linear interpolation operator associated with Q , one has 

ll^'llo,/ < IMIo,/ Vz; € V h (f) , 



where the constant in the inequality depends only on the parameters A and B. 
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Proof. Let Q = : < j < p] and T = {Ij : 1 < j < p}, with hj — \Ij\ = — Cj-i- Let us set 



hi 



if i = 



MG) = < ft* + if 1 < j < P - 1 , 



(6.17) 



h„ 



if j =P 



~ Vi-i- 



Similarly, let Q = {rji : < i < q} and T — {/j : 1 < i < q}, with hi = \Ii 
Furthermore, let h(r]i) be defined in a manner similar to (6.171. Given any v € Vh{T), we have 

imis,/ - E iMfo - E ( w2 fe) + u2 fe-0) ^ = E » 2 &)fc&) . 



and 



Now, for any £j there exist 77, and 8 e [0,1) such that v(£j) = (1 — 6)v(i]i) + 6v(i] i+ i), whence 
w 2 (£j) < v2 (Vi) + v ' 2 ( r li+i)- If # G (0, 1), then £j G (77^,7^+1), hence both and Jj+i intersect 7^. By 
the assumption of locally uniform equivalence of the two grids, we obtain \Ij\ < \Ii\ and < \h\, 

whence h(£j) < hi < min \h{rji), h(r]i + i)^. On the other hand, if 9 — 0, i.e., £j = 77^ , then intersects 
and intersects Ii+\ (with the obvious adaptation if ^ is a boundary point), thus < |Jj| 
and < which yields < h(rji). This completes the proof. □ 

From now on, let I p , lh,p, ^h,D,p, and ICh,D,p, resp., be the interpolation operators defined in 



(|6_5|, fl6.15| and ( |6.16| , resp.. 

Property 6.9. For any p > 0, £/ie operator T p satisfies 

\Ipv\ij < , Vt; g #!(/) , 

tot'/j/i a constant that does not depend on p. 



Proof. The result is classical (see, e.g. 
observing that X p v = T p (Ih, p v). 



). It can be derived from Property 6.4 and Lemma 6.7 



□ 



Property 6.10. Remark 13.5] Assume that cp < q < p for some fixed constant c > 0. Then 

llVllo,/ < IMIo,/ Vt; G P p (J) , 
with a constant depending on the proportionality factor c but not on p. 
Property 6.11. Assume that cp < q < p for some fixed constant c > 0. Then, 

\Zh,qV\ m ,l < \v\ m ,i Vtj G P P (J) , m = 0, 1 , 
tot'/j/i a constant depending on the proportionality factor c but not on p. 

Proof. For m — 0, we observe that Ih, q v = Ih, q {T q v), so that HZ/j^Ho,/ % \\^qv\\oj < IMIo,i by 
Properties |6.4| and |6.10| The result for m = 1 is included in Lemma |6.7| □ 

Property 6.12. Assume that cp < q < p for some fixed constant c > 0. Then, for m = 0, 1 one has 



\ZqV\ m ,I 



\Zh,qV 



qU\m,I 



< 



Mm,/ Vt> G Vh,D,p{I) , 



I^D.^L,/ < Mm,/ Vu G T4, C , p (I) , 
|Zh,C,g«|m,J ^ Mm,/ Vf G Vft, p (J) . 
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Proof. The results for m — follow from Lemma |6.8| applied in various combinations to the grids 
Gq(I), Gp(I), D q (I) and T> P (I), that are locally uniformly equivalent to each other by Corollary 6.6 
The results for m = 1 follow again from Lemma |6.7| □ 

Combining the bounds in the two previous properties, we obtain the following result. 

Property 6.13. Assume that cp < q < p for some fixed constant c > 0. Then one has 

\fCh,D, q v\ m< i < |v| m ,j VuePp(-f), m = Q, 1. 

6.4 Tensor-product interpolation operators 

Tensorization of the univariate interpolation operators considered above yields multivariate inter- 
polation operators on (hyper-)rectangles. Given any such element R = X^ =1 Ik € 1Z as defined in 



fl2.7| ), let 

d 



T P =T« : C%R)^Q P {R), X* = £)2& (6.18) 



fe=i 



be the global polynomial interpolation operator at the nodes of the tensor product LGL-grid G P (R) 



introduced in (2.8 1. On the other hand, let T P (R) = T{G P {R)) be the partition of R into (hyper- 



)rectangles S — Sg = xf = i-ffc,4 for t £ X^_i{l, . . . ,Pk}, where each Ik,i k is an interval in the 
partition T Pk {Ik) of Ik- Then, we set 

d 

V h , p (R) ={»€ C°(R) : v ls G Qi VS* G T P (i?)} - (g) VWI fc ) , (6.19) 
and introduce the piecewise multilinear interpolation operator on the LGL-grid Gp{R) 

d 

l h , p = I*, : C°(R) -> K, P (i?) , l£ p - (g)^ Pfc ■ (6.20) 

fc=i 

Similarly, let T> p (R) = T(Gd,p(R)) be the partition of R into dyadic (hyper-)rectangles E = E m = 
Xfe =1 Dk.m k , where each -Dfe, mfe is a dyadic interval in the partition V Pk (Ik) of Ik- Then, we set 

d 

Vh, D , p {R) = {ve C°(R) : v lE € Qi V£e P p (i?)} = 0Fi, flft (4) , (6.21) 

fc=i 



and introduce the piecewise multilinear interpolation operator on the associated dyadic grid Gd,p(R) 

d 

Ih,D, P = l£ D , p ■■ C°(R) -4 V h , D , p (R) , X^ )P = (g)^ D , Pfc ■ (6.22) 

fe=i 

Finally, we set 

lCh,D,p — K-h,D,p '■ C°(R) -> Vft,x),p(-R) , ^-h,D,p = ' I h,D,p ' I h,p ■ (6.23) 

Given any i-dimensional facet F 6 J|, the analogous definition of the finite dimensional spaces 
Q p (F), V/j^F) and Vh,D,p{F), as well as the corresponding interpolation operators Zjf , lj^ p and 
l£ d p , is straightforward. 

The analysis of our tensor- product interpolation operators relies on the following classical result. 
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Proposition 6.14. For 1 < k < d, let Wk, Vk C be finite dimensional subspaces. Let 

Lk '■ Wk — > Vk be linear operators satisfying, for m = 0, 1, 



l-Mk-r* ;$ \Hm,i k v« e w k . (6.24) 

W := ®fe =1 Wfc an 

/or m = 0, 1, 



Then, setting W :— ®f =1 Wit anrf V :— &>k=i Vk> the operator A — ®f =1 Ak : W — > V satisfies, 



\Lv\\m,R<\Hm,R Vv t W . (6.25) 



/n addition, if the norms in (6.241 /or m — 1 can be replaced by seminorms, i.e., 

\L k v\i Jk < \v\ 1Jk VveW k , l<k<d, (6.26) 



then, the same occurs in (6.25) for m — 1, i.e., 

Lv\ hR < \v\ hR WveW. (6.27) 

The first consequence of this result is the multidimensional version of Property |6.4| It states that 
the operator induces a topological isomorphims between Q P (R) and Vh. p (R) equipped with the 
L 2 or H 1 norms, whose inverse is induced by Ip. 

Property 6.15. For any v £ Q p (R), set Vh :=X^ p (v). Then, 

\\v\\o,R — \\vh\\o,R and \\Vv\\o,r ~ HVw/dlo.fi • (6.28) 

The constants in both relations are independent of p and H. 

Using Lemmas |6 . 7| and |6 . 8 [ Proposition |6.14| yields the following general result. 

Property 6.16. For 1 < k < d, let Qk and Qk be ordered grids in lk, with associated parti- 
tions Tk and Tk, which are locally (A, B) -uniformly equivalent; let Tg k = Tg k : ^{Ik) — > Vh(Tk) 
be the piecewise linear interpolation operator associated with Tk- Consider the spaces Vh(T) ■= 
&>k=i Vh(Tk) and Vh(7~) := ®f=i ^ft(Tfc) of piecewise multi-linear functions on the cartesian parti- 
tions T '■= X^ =1 7fc and T := X^—j Tk of R. Then, the piecewise multilinear interpolation operator 
Zg =Zg ■= ®fe=i^ : C °( R ) V h(T) satisfies 

\\Zgv\\m,R < \\v\\ m ,R Vv € Vh(f) , m = 0, 1 , 

where the constant in the inequality depends only on the parameters A and B. 

In the sequel, the inequality q < p between two multi-indices p, q £ N d is to be understood 
componentwise, i.e., qk < Pk for all k. From Property 6.12 and Proposition 6.14 we immediately 
get the following multidimensional result. 

Property 6.17. Assume that cp < q < p for some fixed constant c > 0. Then, 

\ZqV\m,R ^ \v\m,R Vw <E V h ,D,p{R) , TO = 0, 1 , 

with a constant depending on the proportionality factor c but not on p. 

Finally, from Property |6.13| and Proposition |6.14| we immediately get the following multivariate 
result . 

Property 6.18. Assume that cp < q < p for some fixed constant c > 0. Then, 

\JCh,D,pV\m,R ^ \v\ m ,R Vw € Qp(R) , TO = 0, 1 , 

with a constant depending on the proportionality factor c but not on p. 
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6.5 A localized Jackson estimate for the interpolation error 

For the verification of certain ASM conditions in the next section, we will need Jackson estimates 
for piece wise multilinear interpolation errors, which involve the meshsize in a localized, cell- wise 
manner. Since these results will be applied to both the LGL- and the dyadic tensorized grids, with 
various choices of finite-dimensional function spaces (comprised of either piecewise multi-linear or 
global polynomial functions), we first establish the key estimates in suitable generality in order to 
specialize them later to the cases at hand. 

Consider again the general piecewise multilinear interpolation operator Ig — Ig : C°(R) — > V/,(T) 
introduced in the statement of Property |6.16| above. In addition, assume that each grid Gk is locally 



quasiuniform according to Definition 6.1 For k = 1, . . . , d, let Wk be a finite dimensional subspace of 
to be specified later. Then, the univariate piecewise linear interpolation operator Xg h = Ig* 
is well-defined on Wk and we have 

||Z<fc«||o,j* < CfcNki* Vw e W k , (6.29) 

for some constant Ck > independent of the size |7fc| but possibly depending on the dimension of 
Wk (although this will not be the case in all our applications). Let us set W = (g) fc=1 Wk- 

Next, consider the cells Se forming the partition T = X^ =1 7fe, i.e., T = {Si = X.f,=ilk,e k '■ 
Ik.i k S Tk} , and let hi := maxk\Ik,£ k \ be the largest one-dimensional size of the cell Si. Let 
h = J2i h-eXSe De the meshsize function defined in R. 

The following localized Jackson estimate will be used several times in Sect. [7] for proving certain 
ASM conditions. 

Proposition 6.19. The following estimate holds 

\\h-\v-Igv)\\ , R <\v\ llR VveW, (6.30) 
where the constant implied by the inequality is independent of the meshsize but depends on the 



constants introduced in (6.291 



Proof. The result will be obtained by assembling local estimates in each cell Si, which in turn are 
derived by a scaling argument from corresponding estimates on the unit box B d = [0, l] d , with 
B = [0, 1]. To this end, let I B k := Z-f denote the multilinear interpolation operator on B k , i.e., 

X B uv= w (£) $ ?> 

where denotes the multilinear Lagrange basis function satisfying < &^(^') = for £, £' £ J r (B k ). 
We make heavy use of the fact that I B d is a tensor product operator, i.e., we have I B d — ® d I B . 
Moreover, it will be convenient to employ the following convention to write, for any k — 1, . . . , d, 

B d = A k x B d ~ k , 

meaning that A k is the unit fc-cube representing the first k variables and B d ~ k is the unit (d—k)-cube 
for the coordinates d — k + 1, . . . , d, . 

Lemma 6.20. Let W = Wk, where Wk are finite- dimensional subspaces of H 1 (B). Then, 

one has 

d 

\\v-l B *v\\l Bd <J2 Y, IK^'OIIIU™ . VveW, (6.31) 
fe=i cea^- 1 ) 

(with the first summand on the right-hand side to be understood as \\d Xl v\\^ Bd ), where the constant 
depends only on d. 
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Proof. Notice first that for e := v — I B v one has (I B v)' = v(l) — v(0) — f v'(s)ds so that e(x) — 
Jq v'(s)ds — x f Q v'(s)ds from which one easily derives that 

\\v-I b v\\o,b<2\v\ 1>B) v £ H 1 (B). (6.32) 
Next, denoting by id the identity operator on B, and writing 

d 

v - l Bd v = E <8> id d ~ k+1 )v - (l Bk <g> id d - fe )w) 

fe=i 

d 

= E ( idfe_1 ( id - J s) ® id rf - fc )(I Bfe -i <g> id rf - fc+1 )i;, 
fe=i 

we note that, when abbreviating Wk '■= (X B k-i <£> id d ~ fe+1 u), for x' fe := (xi, , . . . , Xfc-i, £Ek+ij • ■ • , Xd) 
one has 

(id*- 1 <g> (id -J B ) <g> id^JiOfcCsi, . . . , ar fe _i, ifc+i, ■ ..,x d ) = 0, & = 0, 1, 4 e A^ 1 x 

Hence, we can apply ( |6.32| to conclude that 

|| (id*- 1 ® (id-l B )®id d - k )w k \\ Bd < 2||6l e »«'fc||o,B*- 

For k>2, notice that Wk(£, •) = v(£, •), £ £ J r o(A k ~ 1 ) and that Wk is affine in the first k — 1 variables. 
Hence, we conclude that 

B d-k + l , 

from which the assertion of Lemma |6.20| easily follows. □ 
To complete the proof of Proposition |6 . 1 9] consider now a generic cell Si £ T, which we can write 

as 

k d 

Si=(X Iia) x ( X Iia) =■■ A\ x B d f k . 

1=1 l=k+l 

Given any v £ W and considering its restriction to Si, we write v(x) = v{x) with x £ B d , so that 
(Xs e v)(x) = (T B dv)(x); setting hi/, := a standard affine change of variables yields in view 



of (6.311 



\\v-ls e v\\o,s e = \St\P - Zb*Q\\1,b<> ~ I s *! E 51 ll^ fc «(C,-)|lo,B<i-*+i 

fe=i«e^o(^ fc - 1 ) 

< \St\Y, E I^V^JIM^)!!^™ (6.33) 
fc=1 €e^o(A f fc - 1 ) 

<^ 2 E E i^ 1 iii^«(^-)Ho,b™- 

Dividing both sides by hf and summing over i provides 

\\h-\v-i g v)\\l R <Y; E E I^IR^OIIjU™ • 
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Now, \A\ 1 | = Yli=i hi,t t — IIf=i w i-ii f° r eacn £ e -^oO^* L ), where the weights u;;^ are defined 



rfc-l 



fc-l\ 



by the conditions 



E « 2 teKc, = Ifefllo,/, Vi;eC°(/ ; ) 



Then, the assertion follows from (6.291. 



□ 



7 Stage SE-CG -> DFE-CG 

We turn now to the second stage of our iterated auxiliary space preconditioner, namely the construc- 
tion of a preconditioner for the conforming spectral-element problem introduced in Sect. [5] The new 
preconditioner will be based on an additional auxiliary space, comprised of globally continuous 
piecewise affine elements on conformingly matched dyadic meshes. The latter are associated with 
the LGL-grids used in the spectral-element discretization according to the construction mentioned 
in the previous section. 



7.1 Definition of the ASM ingredients 

The roles of the spaces V and V appearing in the ASM (see Sect. K?J) are now as follows. Let us set 



V := V 5 n H^(Q), where V S is defined in (2.151, while V := V KD H H^(fl), where 



V h<D = {ve : VR e TZ , v\ R := v R G V h , D , p (R)} . (7.1) 

Note that since the dyadic grids are by construction nested, the trace spaces obtained by restricting 
Vh.D to any interface of two adjacent elements are non-trivial. Also note that now V <£_ V so that 
the full ASM assumptions have to be checked. 

Since all our finite dimensional spaces are contained in Hq(£1) we set 

a(u, v) = a(u, v) = a(u, v) = (Vu, Vw) ,n = E Vw )o,h = : E a «( u > v ) ■ 

Ben Ren 

We next have to specify the auxiliary form b. A central point of this section is the fact that we 
can no longer choose the form b as in Definition |5.1| In this original form it was designed to absorb 
the jump terms which are now missing. In the current situation it would therefore certainly ensure 



(3.3 1 but would actually be too strong an overestimation of the bilinear form a caused by using 



an inverse estimate on strongly anisotropic hyper-rectangles. As a consequence, one now runs into 



difficulties verifying the conditions (3.6) and ( |3.7| . Our modifications are guided by the following 
closer look at the form a. The restriction a R {u, v) of a(u, v) to any clement R = X Ik is the sum 
of d bilinear forms, i.e., a R (u, v) — '^2'l=i a R,k(u,v), where each form OR,jt is the tensor product of 
univariate forms 



aR,k(u,v) = J d k ud k vdx = J ^ d k ud k vdx k J dx', 

with R' k = X^kh- 

The idea is to "weaken" the bilinear form b by avoiding the inverse estimates in regions where the 



LGL-grid is too anisotropic while preserving the validity of (3.3 1. Our ansatz for the form b(u,v) 
follows the above structure of a, i.e., we set 

d 

b(u,v) = E b R(u,v) , b R (u, v) = ^2b R , k (u,v). (7.2) 

Ren k=i 

Moreover, recall that for v £ Q P (R) and Vh ■= T R p (v) our norm equivalences imply a Ry k(v,v) ~ 
o-R k( v hi v h)- Hence, it suffices to arrange that a, Rt k( v hi v h) S b R .k(v,v). We proceed now defining 
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fr-R.fcC^i v) in a cell-wise manner for any u, v £ V + V. To that end, we first decompose R into the 
ci-dimensional LGL-subcells Se — Sg(R) = X fc=1 ifc^ for ^ G X f = i{l, ■ ■ • ,Pfc}, where p — p{R) and 
each partition T k (h) = T{G Pk {h)) = : 1 < * < Pk} is generated by the LGL-grid G Pk {h)- 

Next, for each 1 < k < d, we decompose the partition T{R) = {Se(R)} of R into two parts T k °\R), 
T^(R) defined as follows. Setting hi t t t := \ and fixing a constant C aspoc t > 0, we let 

S e eT£ 0) (R) if ^y^'^ > aspect, S e eT, (1) (i?) otherwise. (7.3) 

Thus T k °\R) and T^\R) are comprised of "strongly anisotropic" and "sufficiently isotropic" cells, 
respectively. In analogy to a R ^ k (u h , v h ) = X^eT(.R) a R,k,s e (u h , v h ) we make an ansatz for b R , k (u, v) 
as 

6*, fc (u,«) = X! CaK^ E b R,k,s e ( u ^) =■■ bf, k (u,v)+b%(u,v), (7.4) 

so that we still have an t k,s e ( u h,Vh) ^ g (u, v) for all G 7% (R), i = 0,1, as follows. In 
analogy to the factorization R = X R' k , we define for every Se G 7~(R) the (d — l)-dimensional 
hyper-rectangle S' gk := S' t k (R) :— X J=li ;^ fe Ii t i v The idea is to bound the terms 

d Xk u h d Xk v h dx = I ( / d Xk u h (x k ,x')d Xk v h (x k ,x') dx k ) dx' 

Si 

by using quadrature in all but the fc-th variable while keeping the integral with respect to the k-th 
variable for Se G T k °\R), and applying an inverse estimate for Si G Ti , (i?). Specifically, using the 
tensorized trapezoidal rule for integration over S' g k , which is the finite-element lumped mass matrix 

approximation 15, (4.4.44) on p. 220], we define for St G T fc (0) (i?) 

b< R,k,s e ( u ' v ^ := E w « / d Xk u h (x k ,C)d Xk v h (x k ,^)dx k , (7.5) 

where the weight 



d 

J 'e. k 



II hi,e t — vo ^-d-i(Si k ) (7.6) 

is equal to the (d— l)-dimensional volume of S' e k , and as before, FoiS'g k ) is the set of all its vertices. 
Note that this set is isomorphic to G p (r)(R) H (S' e k x {£}), where £ denotes any cndpoint of the 
interval J Mfc = [6i,4,6,4+i]' 



On the other hand, for each Se G Ti (R) we apply to each summand of (7.5 1 the inverse estimate 



[d Xk u h (x k ,x')] dx k < — [«fc(6Mk> x + u fc(&A+i> a; ')] » 

which leads us to define for St G Tu^ (R) 

C t ( tt '") := E E ^'^-«,,(f,o«fc(e.o. a7) 

£'e^o(sj <fc )££^o(iiM fc ) Mfc 

where as before the c$ £/ are tuning constants of order one which have to be chosen judiceously in 
practical applications. To simplify the exposition we shall suppress them in what follows. 
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In summary, we have 

&(«,«) = ££( £ &£U («.«)+ £ (7 - 



where K'. s (u, i = 0,1, are defined by (7.5 1, (7.7), respectively. By construction this bilinear 



form satisfies the ASM condition (3.3 1 
Property 7.1. One has 



a(v,v) < b(v,v) VveV. (7.9) 



Note that the relation (7.9 1 can be viewed as an inverse estimate. 



Let us define next the operators Q and Q, starting with the former one. Obviously, this can be done 
element-wise, so let us fix again an element R G 1Z. Given any v G V, let us set vr := v\r G Vh,D,p{R)- 
Let us denote for any vertex z G of R by 

{E 1 ,E 2 ,...,E d }=F 1 {z)nF 1 {R) 

the set of edges of R containing z, with E^ parallel to the fc-th coordinate direction. Let us also 
introduce the vector of polynomial degrees 

P * z = (j>t(E 1 ), p t(E a ),...,pl(E d ))eN d 



and note that, by definition of p^(Ek) (recall ( |5.8| )), we have p* z < p{R) componentwise. Finally, 
let us introduce the localizing function Q z G Q 1 (R), defined by the conditions & z (y) = 8 VtZ for all 
y G J-q(R). Then, we set 

K-=^h,D, p; (^zV R )eV ht n, Pt (R) and < = 2* v* € Q p * (R) . (7.10) 

Summing-up over the vertices of R, we define 



V R 



£ w* G Vh,D,p{R) and Q R v R :=v* R := £ v* G Q P (R) ■ (7.11) 
z£F (R) z<E^ (R) 



The following properties will be useful in the sequel. 
Property 7.2. For any edge E G J-\(R), one has 

(v*r)\e = (vr)\e ■ 

Proof. If E has vertices z\ and z 2 , then 

+*U)\E =I h,D,p»(E) + |JS =Ih,D,pt(E) (*r)\E = (««) \E 

since by continuity (vr)\e G V^ jj ^^e^E); on the other hand, (v*)\e — for all y G J-q(R) \ {z±, z 2 }. 
Thus, = (v* Zi +v* Z2 )\ E = (vr)\e- □ 

Property 7.3. For any interface F G T&—\, with TZ(F) = {R',R"}, one has 

{vr>)\f = {v*r»)\f and {v* R ,)\ F = {v* w ,)\ F . 
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Proof. Let J-o(F) = J-o(R') n J-q(R") be the set of vertices of F. Since, by continuity of v, vp := 
(vr')\f = {vr»)\f, one has ($ z v R >)\ F = ®gV F = ($ z v R »)\ F , z G Fq(F), where $f G Qi(F) satisfies 
&z(y) — Sy^ z for all y G Fq(F). Hence, denoting by q* z — (p*)' G N^ 1 the reduced degree vector 
obtained from p* by dropping the component in the direction orthogonal to F, one concludes that 



F v F ) 



zer (F) ze^oC-F) 

= E X h,D,qi {® Z VR»)\F = {V*R»)\F , 
zeF (F) 

where we have used that q* is the same for both R' and R" . This confirms the first part of the 
assertion. Abbreviating vf := Tf lD » [^vf), the second one follows from 

i v R')\F = E X 9i = ( V *R")\F , 
z€Fo(F) 



which completes the proof. 



□ 



The previous property guarantees interelement continuity of the functions v* R and v* R , justifying 
the following definition of the operator Q. 



Definition 7.4. Let Q : V — > V be given by 

(Qv)\r = Qrvr = v* R 



V-R G TZ 



where v* R is defined in (7.11) 



The definition of the operator Q follows the same lines as above. Given any v e V and any R € TZ, 
we set vr = v\ R G Q P (R). Then, the chain ( 7.10| -( |7.11| is replaced by the following one: 



v* :=1* eQ Pl (R) and 

Summing over the vertices of R, we define 



v* R := E v l e Qp( R ) and Qrvr '=Vr= E v* z S K,_D, P (-R) 



(7.12) 



(7.13) 



As above, one easily confirms interelement continuity, which suggests the following definition. 
Definition 7.5. The operator Q : V — > V is given by 

(Qv)\r:=QrVr = v* r VR G TZ Vu G V , 



where v* R is defined in (7.13). 



7.2 Check of the ASM assumptions for Q and Q 



Hereafter, we deal with the ASM conditions (3.5 1 and (|3.7[) involving the operator Q, as well as (3.4 1 



and (3.6) involving the operator Q 



Proposition 7.6. The operator Q is linear and satisfies the continuity assumption (3.5 I 
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Proof. We have to prove that |Qt>|i,n ^ Hi,fi for all v £ V. This follows if, for any R £ TZ, we prove 
that 

K\i,R < \vr\i,r Vv r g V h . D . p {R) . (7.14) 
A classical mapping-and-scaling argument in Finite Element analysis tells us that this result holds 
provided it holds when R is the reference element R = Xj_x /, with I = [—1, 1]. For this element, 
it is enough to prove that 

\ v \r^\Hi,r VveV hiDtP {R). (7.15) 
Indeed, changing v into v + c, by any ceM, does not change the left-hand side, whence 



In order to establish (7.15), let us first consider a univariate function v £ Vh,D,p{I) an( i let $ denote 
the affinc function taking the value 1 at one endpoint of the interval and at the other one. Let us 
prove that if cp < q < p for some fixed constant c > 0, one has 

Ph,D, q (*Z)\ m ,t<\M m ,t, m = 0,l, (7-16) 
where, of course, the involved constant depends on the proportionality factor c. For m = 0, we have 

\\ik,D, q m\\h < £ mcH0fhn, q (0< £ %) a w»?) s \\ni, t > 

since <I> 2 < 1 and hn,q(0 ^ ho.piC) for all ( £ Go,q(I) C Gd,p(I)- For m = 1, we have by Lemma 

\ih,DA^)\ij^\^\u^Mi,v 



6.7 



where the last inequality holds since we are working on the reference element. Hence, (7.16) is 
proven. Using this result and Proposition |6.14| we obtain the bound 



for each function v* defined as in (7.101 on R. Then, Property 6.17 yields \vl\ 1 R < 
(7.15) follows by the triangle inequality, since |J-o(-R)| — 1. 



and 

□ 



Proposition 7.7. The operator Q satisfies the Jackson condition (3.7l ; i.e., one has 



b(v -Qv,v- Qv) < \v\f a 



Wv£V 



where the multiplicative constant in the above estimate depends on the constant C aspcc t in (7.3 I. 



Proof. For each R £ TZ and each k = 1, . . . , d, let us recall the definitions ( 7.4 1 , (7.5| and (7.7! of the 
form bR^{u,v) and its portions b^ k (u,v) and b^ k {u,v). Concerning the former portion, observe 
that 



b Rk( v > v )= £ £ w 'i,k I \d Xk Ih,p v (C,x k )\ 2 dx k 

£ I \d Xk l*: p v(x)\ 2 dx < £ I \d*A<x)\ 2 dx<K P v\iR, 

.fn\ . J Sf> c •~r~ { n i <J Sf 



(7.17) 



s e eT k m (R) 



S e eT(R) 



where we have used the uniform equivalence of the weights uo' lk (see (7.6 1 ) with the integration 
weights for multi-linear functions on the cell S' e k . Thus, for all vr £ Vh.D.p{R), we have 



b R,k(VR - QrVr, Vr - QrVr) < \T^ p {v R - QrVr)\\r 

< K p {vr)\Ir + K p {Qrv R )\Ir < \vr\Ir 



(7.18) 
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where the last bound follows immediately from Lemma |6.8| Properties |6.11| and |6.16| and Proposi- 
tion ESI 

Consider now the portion b R \(u, v). Arguing as above, one has 

O v ' v )= E E < k E ^Ii^X^oi 2 ^ E / Ki\i« p v(x)\?d X . 



S e £T k {1) (R) 



Now observe that, by definition of 7^ (i2) (see (7.3 1), one has h k ^ k > max(l, C aS p ect ) hi, with 
hi = max; hi g, . Hence, 



£ f hj*\l* p v(x)\ 2 dx< ]T [ hf\ll p v{x)\' 
s e er k w (R) Se S e €T(R) JS * 

,-1-T-R „,l|2 <- i-R „,M|2 , IU-L.II2 



= \\h- L I^v\\i iR < ||/i- 1 («-^ p «)||5,H-r||i. «iio,h , 
where h = h^xse is the LGL mcshsizc function in R. If G V/ l; £) jP (i?), this yields 

b R , k (v R - Q R v R ,v R - Q R v R ) < \\h^ 1 (v R -2^ p v R )\\l R + \\h~ 1 (Q R v R -l^ p {Q R v R ))\\l < 

+ \\h~ 1 (v R - QrVr)\\q R ■ 



(7.19) 



Now we invoke Proposition |6.19| with different choices of the grid Q and the space W, to bound 
each of the three summands on the right-hand side. For the first summand, we have Q = G P (R) (the 
LGL-grid of order p in R) and W — V} ltR>tP (R). We note that, due to Lemma 6.8 the bounds (6.291 
are satisfied with c& < 1. Thus we get 

\\h-\v R -X^ p v R )\\l R <\v R \l R . (7.20) 

For the second summand, we recall that Q R v R G Q P (R), so that we have again Q = Q p (R), whereas 

" Thus, 



now W — Q P (R ). The bounds (6.291 are now satisfied with < 1 because of Property 
recalling (7.141, we obtain 

\\h-\Q R v R -1% p (Qrv r ))\\I r < \Q r v r \Ir < \MI,r ■ 



6.11 



(7.21) 



At last, we bound the third summand on the right-hand side of (7.191. To this end, recalling the 
definition (7.111 of Q R v, noting that v R = X) z eJ r (fl) ^zV R , and defining v z := $ z v R , z G ^(R), it 
suffices to bound the quantity ||/i _1 ('Dz — Zp* (%hD P *^z))llo R f° r eacn z e J~ (R). Writing 



and setting for simplicity w z = £> p*^ e ^ft,,r>.p* we have 



R „7, 



(7.22) 



We now proceed as in the proof of Proposition 



7.6 



i.e., we work on the reference clement R viewed 



as an affine image of R. This simplifies handling the factor <& z when eventually bounding the H 1 - 
scminorm of v z — <fr z v R by that of v R on the element R. 

We want to apply Proposition 6.19 to the first summand on the right-hand side of (7.22), with 
Q = V p *(R) (the dyadic grid of order p* in R) and W = & z Vh y D, P *(R)- To this end, we o bser ve 
that the meshsize function h associated with the grid G P (R), as a consequence of Corollary 



6.6 



uniformly comparable to the mcshsizc function ho, P ' associated with the grid T> p * (R), i.e., h ~ ho, P * 
pointwise in R. On the other hand, the bounds (6.291 are satisfied with c k < 1; this easily follows 
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from the fact that the restriction of v z to any cell of the grid T> p * (i?) is a piecewise multi-quadratic 
function belonging to a finite dimensional space whose dimension is bounded independently of p. 
Thus, we obtain 



\h-\v s -I% D j.v z ))\\ 



. < 



0,R 



J R"l,R ■ 



(7.23) 



The second summand on the right-hand side of (7.22 1 can be bounded with the aid of Proposition 6.19 



as well. Indeed, using the property that Xh,q(Xqv) — Xh. q v \iXh q v and X q v are the low- and high-order 
interpolants of a continuous function on the same grid, one has 

v - X q V = V - I h , q (XqV) + X Kq (X q v) - X q v = (v - X h , q v) - (X q v - X h ^ q (X q v)) . 
In our situation, this yields with u z — X^,w z £ Q P *(R) 



\\h-'(w z -X*w z X A < \\h-'(w z -X^ p ,w z W oA +\\h-'(u z -X^ p , 



jR f, Ml2 



0,R 



so that we can apply Proposition 6.19 with Q = G P *{R) and either W = Vh,D, P {R) or W = Q P *(R). 
Again, the function h is uniformly comparable to the function h p * associated with the grid G P *(R), 
and one easily checks that the bounds (6.29 1 are satisfied with Cfe < 1 with both choices of W. Thus, 



\\h-\w z -X*w z )t Qfi < \w z \l A < \v z \l A < \\vJl 



(7.24) 



where the second inequality follows from Lemma |6.8| and Property |6.16| Going back to the element 
R, the bounds (7.221), (17.231) and \7.2ty imply 



\\h-\v z -X^{X^ p ,v z ))\\l R <\v R \l 



R ■ 



which, together with (7.20) and (7.211, allows us to obtain from (7.191 

b R,k(vR ~ QrVr,v r - Qrv r ) < \v R \\ R . 



This completes the proof of Proposition 7.7 



□ 



In a similar manner as above, one can check the ASM assumptions for Q, where essentially the 
order of the operators X p ,Xh,D,p is interchanged. However, Lemma 6.8 allows us to argue as before. 



Proposition 7.8. The operator Q satisfies the continuity assumption (3.4 1 and the Jackson as 



sumption (3.6 1 



So far we have shown that the above choice of auxiliary spaces complies with the ASM require- 
ments so that, in principle, the stage II preconditioner is again optimal. However, again it remains 
to precondition the new auxiliary problem, i.e., to identify C^. As the auxiliary space V now 
corresponds to a hierarchy of nested low order finite element spaces, this task looks more feasible. 
Nevertheless, the strong anisotropy appearing in the underlying dyadic meshes poses a further ob- 
struction. In fact, an asymptotically optimal treatment that works for arbitrary degrees seems to 
preclude the application of standard tools. Therefore, we postpone this issue to the forthcoming 
part II of this work where a multi-wavelet preconditioner will be proposed and analyzed. Here we 
address in the following section a second obstruction, namely the fact that the Gramian B associated 



with the auxiliary form b, defined by (7.8 1, is no longer diagonal, due to the coupling across clement 



interfaces and the remaining integrals of derivatives in strongly anisotropic grid areas. 
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8 Numerical results 

In this section we quantify the performance of the multi-stage preconditioner, examining first the 
influence of the various ingredients in each of the two stages discussed in this paper. This includes the 
resulting condition numbers and their dependence on the parameters, in particular on the polynomial 
degrees. In the experiments the condition numbers are obtained as quotient of the largest and 
smallest eigenvalue of the preconditioned linear systems of equations. Since the matrices involved 
are large and sparse, we apply a Jacobi-Davidson type method for the eigenvalue computation that 
is tailored to the specific needs in this application, see [5] for details. 



8.1 Test cases 

We now consider two test scenarios that represent typical situations involving varying polynomial 



degrees when solving a Poisson problem (2.13) on the domain fl = [0, 3] 2 C M 2 that is divided into 



nine square patches of equal size as depicted in Fig. 2(a) This grid of 3 x 3 patches contains the 
central patch that is not in contact with the domain boundary. Inside each patch we use at this point 
just constant degree vectors, i.e., the polynomial degrees are the same in x- and y-direction inside 
each patch. However, the polynomial degree may vary from patch to patch and we are interested 
in quantifying the influence of patchwise varying polynomial degrees on the performance of our 
preconditioning technique. 

In this first test scenario we arrange the polynomial degrees in the nine patches to be either p or 
q in a checkerboard-style pattern. Concerning the values of p, q G N, we consider the following five 
cases: (i) we either use constant polynomial degree q — p, (ii) simulate a small variation of the degree 
by choosing q = p + 2, or a large variation of the degree represented by multiplicative relations (iii) 
q = 3/2p for even p, (iv) q = 7/4p for p chosen as a multiple of 4 or (v) q = 2p. 
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(a) First test scenario: checkerboard-style grid. 

Figure 2: Grid and distribution of polynomial degi 
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(b) Second test scenario. 



for the test scenarios. Darker shading indicates 



In the second scenario we simulate a situation that typically arises in p-adaptation towards a 
singularity, namely we solve Poisson's equation (2.131 on f2 again split into 3x3 square patches. 
Assuming a singularity located at the origin, the polynomial degree increases in directions away from 
the singularity as illustrated in Fig. |2(b)| 
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8.2 Stage SE-DG -> SE-CG 

In the first stage, we apply the auxiliary space method to precondition the DG-spectral element 
problem by an auxiliary problem given by a conforming Galerkin formulation. We first investigate 
the resulting condition numbers and their dependence on the polynomial degrees for that stage only, 
i.e., the auxiliary problem is solved exactly. 



8.2.1 Choice of the bilinear form &(•,•) 



In Subsect. |5.2| the auxiliary bilinear form b in the first ASM stage is defined up to "tuning" constants 
C£ ~ 1 that still need to be specified for £ g G P (R), R E 7Z. In view of the inverse estimate used in 



the proof of Proposition 5.6 we define b$ : Vs x Vg — > K as 



^m=akE E «(0«(0w' e +7p 1 E w ^E E . 

where Wpr± is the LGL quadrature weight on F seen as face of R ± . Consequently, a reasonable 



ansatz for the constants in (5.7 1 is 



_ / Pi(4 + lPiu F w F , R /W(.), for £ G G P (F, R), F g 7" d _i(ii), R 
PiCi, else. 



Preliminary experiments concerning the constant arising in the inverse estimate reveal that c\ = 10 
is a good choice which we fix in our subsequent tests. 

We emphasize that with this choice of b the corresponding matrix B is diagonal for any choice of 
/?i > and p\ > 0, i.e., the complexity of solving a linear system with the matrix B is proportional 
to the number of unknowns. 



8.2.2 Optimal choice of the parameters /3\ and p\ 

We next address the question of suitably specifying the parameters /3i and p\ in the bilinear form 
b. In particular, is there a good choice for these parameters that works by and large independently 
of other parameters, such as the polynomial degree? We investigate this problem for the first test 
scenario and and the (extreme) case (v) q = 2p. The condition numbers kasm = k(CaA) of the 
linear system of equations preconditioned by the auxiliary space method, where the auxiliary problem 
is solved exactly are depicted for p = 8 and for p = 16 in Fig. [3] as functions of (3\ g [0.05, 1.1] and 
Pi G [0,2]; contour plots, i.e. the lines of constant condition number, are displayed. 

We observe that the simple choice fii = 1 and p\ = 0, which is eg = 1, is not optimal. In fact, for 
both p = 8 and p = 16 there is a very flat minimum of the condition number that is located near 
the parameter point (/3i,pi) = (0.15, 1.25). From now on we fix these parameter values for the rest 
of the paper. 



8.2.3 Dependence on the polynomial degrees 

Now we are ready to investigate the key issue, namely the robustness of the preconditioner with 
respect to the polynomial degrees. Fig. |4(a)| shows the condition numbers obtained in the first 
scenario for the relations (i) - (v) between p and q. We observe that the condition numbers stay 
uniformly bounded as p increases although the upper bound depends on the ratio q/p. In the case 
of additively increasing the polynomial degree q — p + 2, the quotient q/p decreases which is also 
visible in Fig. |4(a)| 

The analogous plot for the second test scenario, representing a typical p-adaptation, is depicted in 
Fig. |4(b)] In the second test scenario the condition numbers are almost constant and slightly smaller 
than 7.5. 
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(a) p = 8, q = 16 (b) p = 16, q = 32 

Figure 3: Contour plots illustrating the dependence of the condition number kasm = k(CaA) on 
Px > and Pl > 0. 

8.3 Stage SE-CG ->■ DFE-CG 

Now we analyze the performance of the preconditioner in the second stage. 
8.3.1 Solution of the smoothing problem 

Every application of the auxiliary space preconditioner in the second stage requires solving two linear 
systems of equations: the auxiliary problem involving the matrix A and the problem involving the 
matrix B, which we call the smoothing problem. Unfortunately, the corresponding matrix B = B 2 
in the second stage is in general no longer diagonal. The modifications on the strongly anisotropic 
cells in T^°\R), R S 1Z, cause a slightly stronger coupling than for a purely diagonal structure 
although B2 remains diagonal in the bulk of the element. We therefore need an efficient solution 
for the smoothing problem as well. We first consider the case where the domain consists of a single 
patch since this already offers valuable insight into various effects. 



Single patch: sparsity pattern and reordering Regarding the possible use of both (blockwise) 
direct or iterative solvers, we first study the sparsity pattern of B2 obtained when the polynomial 
Lagrange nodal basis is ordered in a line-wise fashion. In particular, since the Lagrange nodal basis 
functions with reference nodes at the boundary of the patch need a special treatment anyway, we are 
interested in the sparsity pattern of the submatrix for ansatz functions in the bulk of the patch. In 
Fig. 5(a) this is depicted for a patch with polynomial degree p = 25 in both directions. We confine 
the discussion in this paper to the bivariate case d = 2. 

Since the cells of high anisotropy are selected on a per-subcell basis the inverse estimate is used 
wherever possible. Therefore, B2 is still very sparse. The two bilinear forms s and bp^ 2 s , 
see (7.5 1 and (7.7) for their definitions, contribute only to the main diagonal of the matrix. Since 
a rectangular cell in two spatial dimensions can only be anisotropic in at most one direction for 
Caspcct > 1 the intersection of the regions where s and Si are applied does not contain 
a subcell. At most this intersection might be comprised of isolated nodes. Because of the quasi- 
uniformity of LGL-grids the regions Tx°\R) and T 2 (R) do not overlap if C a s P oct is large enough, 
see e.g. |8J|5]. 
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Figure 4: Condition numbers kasm = k(CaA) obtained in stage SE-DG — > SE-CG. 
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Figure 5: Sparsity patterns for patch-inner nodal functions: (a)-(b) top left quarter of the matrix 
B2 of p — 25 on a single patch, (c)-(d) Matrix B2 on two patches with polynomial degrees 
6 x 6 and 12 x 12. 



Moreover, due to the mass lumping approximation involved in the definition of the bilinear form 



6, see (7.5) and (7.7 1, one can check that, given two tensorized Lagrange polynomials <f>^ and 



associated to LGL nodes in Q p (R), one has 

M*£>{) + =► 3 fc, 5, G T fc (0) such that £ k ,£ k e Mhi k ) and = | fc , for k' + k. (8.1) 

In other words, there are at most three diagonals populated with nonzero entries in each of the 
matrices corresponding to and \ , but the positions of the diagonals besides the main diagonal 
differ in both matrices. 

In order to bring the non-trivial diagonals as closely together as possible we use a special ordering 
exploiting the fact that for C asp ect large enough the sets of nodes contained in T x (R) and Tz°\R), 
respectively, are disjoint and can each be reordered separately. More precisely, we deviate from the 
default line- wise ordering in the first coordinate direction only for patch-inner nodes in 71 (R), 
k 7^ 1, switching to a line- wise ordering in direction k, see Fig. [6] for the details. Note that this 
reordering technique can still be applied for smaller C aspoc t when there are nodal basis functions 
where the corresponding node is contained in 7j (R) n 7~2 (R) as long as these basis functions do 
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not produce an off-diagonal entry with other patch-inner basis functions in 7~i (R) U 7^ (R), see 
e.g. nodes 16 and 19 in Fig. |6(b)| 

The resulting linear system of equations has a block-diagonal matrix of tridiagonal submatrices, 
sec Fig. [5(b)] and can be inverted in O(N) operations, where N is the size of B2. 



(so) (87) (ss) (m) (90) (91) (92) (93) (94) (95) (9g) (97) 

(72) (73) (74) (75) (76) (7?) (78) (79) (so) (s?) (S2) (83) 

(ss) (59) (go) (gi) (02) (oa) (ba) (bs) (gg) (g?) (gs) (go) 

(44) (45) ^y(^y(^y(^)(52^ (53) (54) (55) 

(30) (3?) (32)-(33) (34) (35) (36) ©)(©""©) (40) (4?) 
(l6)-(^(is) (19) (20) (2l) (22) (23) (24) (25) (26)-(27) 

(a) Line-wise numbering. 



(so) (87) (ss) (89) (So) (9?) (92) (93) (94) (95) (9G) (97) 

© © © © © © © © © €) (£) © © © 

(58) (59) (go) (g?) (G2) (g?) (04) (gs) (gg) (g7) (gs) (W) 

© © © © ©j5H©©>©l© © ©©© 

(20) (2?) (23)-(2?) (^9) (44) (47) (50) {©-©) ©1 ©) 

© ©*"©)© © © © © © © ©)©}©) © 
(b) Reordered numbering. 



Figure 6: Schematic illustration of the numbering of the inner degrees of freedom in a single patch R 
of polynomial degrees 13 x 13 for C asp ect = 1.6. Because of the symmetry, only the lower 
half of the patch is shown. Black nodes: boundary degrees of freedom, gray nodes: inner 
degrees of freedom corresponding to vertices of cells in 7^ (R)> white nodes: inner degrees 
of freedom corresponding to vertices that are not vertices of cells in 7^ (R). 



Multiple patches: Gaufi-Seidel relaxation on the skeleton When 1Z consists of several 
patches, the situation is more complicated. Note that although b is defined on V, it is only evaluated 
for piecewise polynomial functions from the space V. Due to the continuity of functions in V 
across cell interfaces, the matrix B2 in general is no longer tridiagonal as shown by Figures |5(c) 
and |5(d)| There the sparsity pattern of the matrix B2 is depicted for the situation of two patches 
of polynomial degrees 6x6 and 12 x 12 which meet at a common interface. The top left part 
of the matrix corresponds to the degrees of freedom of the patch with lower polynomial degree, 
while the degrees of freedom of the patch of higher polynomial degree are represented by the lower 
right part. Since the trial functions are polynomials on each clement R, the coupling is based on 
polynomial interpolation causing the additional non-zero entries off the central bands. When the 
polynomial degrees on both sides of an interface do not agree, every basis function on the element 
with lower degree corresponding to a node on the interface has to be expanded as a linear combination 
of higher degree Lagrange basis functions on that interface. The continuity requirement therefore 
causes that a lower-degree Lagrange polynomial is coupled with a linear combination of face-inner 
Lagrange polynomials on the higher-degree side. At a common element vertex, the coupling involves 
all Lagrange polynomials with reference node located on an interface containing the element vertex. 

In order to solve the linear system with B2, we alternate a Gaufi-Seidel relaxation over the skeleton 
of 1Z and a block elimination for each R £ 7Z in the spirit of substructuring methods. For the 
patchwise block elimination process for the bulk degrees of freedom we invoke the method detailed 
in the preceding paragraph. 



8.3.2 Dependence on the polynomial degree 

As in Sect. |8.2.3| for the first stage, we now analyze the performance of the preconditioner for the 
second stage. The dependence of the condition numbers of the preconditioned linear systems on the 
polynomial degrees is depicted in Fig.[7]for the two test cases. In all cases we choose a — 1.2, which is 
a reasonable choice that balances two effects: on the one hand, the auxiliary space V — Vh,,r>ni?o(f2) 
is rich enough for a good approximation of elements of V, on the other hand, the dimension of V 
is not too high such that the linear systems of equations for the auxiliary problem is not infeasibly 
large. 



38 K. Brix. M. Campos Pinto, C. Canuto and W. Dahmen 



Note that the nodes of an LGL-grid gradually move with growing polynomial degree. In contrast, 
due to their more rigid structure, the associated dyadic grids evolve more abruptly with occasional 
stagnation. This seems to cause some oscillations in constants related to the interplay between the 
LGL- and the associated dyadic grids, see also [TU] , 

Furthermore we have to fix the constant C aspcct that limits the aspect ratios of subcells to which 
an inverse estimate is applied in the formation of the auxiliary form b. In the recorded numerical 
experiments we set C asp oct = 2. For larger C asp ect the condition numbers up to polynomial degree 
p = 40 exhibit a more pronounced oscillatory behavior. Although they appear to lower again after 
a burst of such oscillations a clear asymptotic bound cannot be observed in this range yet. As in 



the first stage in Subsect. 8.2.2 the tuning constant ~ 1 is still a free parameter. Due to the 



oscillations in the condition numbers, it is less clear how to find an optimal choice because we are 
lacking a clear objective functional. In preliminary studies we found that c^i — 0.6 is a good choice 
which partially compensates for the effect of the constant appearing in the inverse estimate. 
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Figure 7: Condition numbers kasm = k(CaA) obtained in stage SE-CG — > DFE-CG 



Figures 7(a) and |7(b)| show the condition numbers obtained in the second stage by the auxiliary 
space method in the first and second test scenario, respectively. The marked points indicate the 
condition numbers obtained by approximately solving the smoothing problem with the aid of 7 
iterations of the substructuring method. In contrast, in both plots the continuous graphs represent 
the results when the smoothing problem is solved exactly using a direct method. 

For the first test scenario the condition numbers are presented for p ranging from 4 to 40. In the 
cases (i) q = p, (ii) q = p + 2 and (iii) q — 3/2p the condition numbers vary only mildly, probably 
due to the discrete process of dyadic grid generation. In contrast, we observe an increase of the 
condition numbers when larger jumps between the polynomial degrees on adjacent elements are 
considered according to the cases (iv) q = 7/4p and, in particular, to (v) q = 2p. A more thorough 
analysis reveals that this phenomenon is caused by the slow convergence to the asymptotic limit 
of the constant in the estimate analogous to ( |7.16 1 for m = 1, where the roles of the LGL- and 
dyadic grids are interchanged. This estimate enters the proof of the properties of the operator Q in 



Proposition 7.8 and affects the dependence of the constants in the ASM conditions (3.4 1 and (3.6 1 on 



the polynomial degree. This in turn causes some growth of the condition numbers for the moderate 
polynomial degrees under investigation. For a detailed discussion of these approximation estimates 
and the convergence behavior of their constants we refer to |10j . 

One observes enhanced oscillations when specific polynomial degrees are present in the grid, e.g. 
for odd polynomial degrees close to 27 which we attribute to a resonance effect between the LGL- and 
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dyadic grids. More detailed investigations reveal that similar oscillations are also observed for some 
even polynomial degrees, e.g. 50. Moreover, the resonance can be shifted to nearby polynomial 
degrees when the parameter a, controlling the associated dyadic grids, is varied. Therefore, our 
tests should be viewed as a guide for favorable choices of the dyadic grid generation parameters a 
depending on the employed polynomial degrees. 

For the second test scenario, which is perhaps more relevant for practical applications than the 
first one, we investigate p in the range of 4 to 60. When the smoothing problem is solved exactly, we 
observe that, aside from some oscillations due to the grid resonance effect, the condition numbers 
grow in essence mildly for small polynomial degrees and quickly tend to a limit below 10. When 
solving the smoothing problem only approximately by a fixed number of substructuring iterations, 
we observe increasing condition numbers for polynomial degrees larger than 41. This issue will be 
examined more closely in forthcoming work. 

For instance, using &(•, -)-orthonormal polynomials on each patch R further sparsifies B2 signifi- 
cantly. Depending on the range of p, the cost of changing bases may well be offset by the positive 
effect on the substructuring iteration. Detailed studies of this and related issues are deferred to 
forthcoming work. 

Of course, from a practical point of view polynomial degrees exceeding 40 significantly are rarely 
encountered except in pure spectral approximations on a single element when the solution is very 
regular, and in this case the auxiliary problem can be solved efficiently with the method presented 
above. 



8.4 Combined stages SE-DG -> SE-CG and SE-CG -> DFE-CG 

After studying the effects of the two stages and their numerical effects in detail, we now combine the 
stages SE-DG — > SE-CG and SE-CG — > DFE-CG by using the second stage as a preconditioner 
C^ of the auxiliary problem of the first stage. The numerical results for both test scenarios showing 
the condition numbers obtained when the combined preconditioner is applied are depicted in Fig. [8] 
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Figure 8: Condition numbers kasm = k(CaA) obtained in combined stages SE-DG 
and SE-CG -4 DFE-CG. 
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For both scenarios, the numerical effects observed for the single stages are also present when the 
combined preconditioner is used. The oscillations that have been striking in the second stage for 
the first scenario now only have a mild effect, but their dependence on the ratio q/p is still clearly 
visible. In the second scenario, for p < 40 the condition numbers for the approximate and exact 
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solutions of the smoothing problem almost agree. For larger p the insufficiently accurate solutions 
of the smoothing problem by the substructuring method in the second stage lead to a deterioration 
of the overall condition number. When the smoothing problem in the second stage is solved exactly, 
the condition number is bounded by 17. 
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